Skip to content

Commit

Permalink
improve triangulation and LearnerND typing
Browse files Browse the repository at this point in the history
  • Loading branch information
basnijholt committed Dec 15, 2019
1 parent f21f19d commit e959e8c
Show file tree
Hide file tree
Showing 2 changed files with 51 additions and 44 deletions.
2 changes: 1 addition & 1 deletion adaptive/learner/learnerND.py
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ def volume(simplex: List[Tuple[float, float]], ys: None = None,) -> float:
return vol


def orientation(simplex):
def orientation(simplex: np.ndarray):
matrix = np.subtract(simplex[:-1], simplex[-1])
# See https://www.jstor.org/stable/2315353
sign, _logdet = np.linalg.slogdet(matrix)
Expand Down
93 changes: 50 additions & 43 deletions adaptive/learner/triangulation.py
Original file line number Diff line number Diff line change
@@ -1,14 +1,18 @@
import collections.abc
import math
from collections import Counter
from collections.abc import Iterable, Sized
from itertools import chain, combinations
from math import factorial
from typing import Any, Iterator, List, Optional, Sequence, Set, Tuple, Union
from typing import Any, Iterable, Iterator, List, Optional, Sequence, Set, Tuple, Union

import numpy as np
import scipy.spatial

Simplex = Tuple[int, ...] # XXX: check if this is correct
SimplexPoints = Union[
List[Tuple[float, ...]], np.ndarray
] # XXX: check if this is correct
Simplex = Tuple[int, ...]
Point = Union[Tuple[float, ...], np.ndarray] # XXX: check if this is correct


def fast_norm(v: Union[Tuple[float, ...], np.ndarray]) -> float:
Expand All @@ -21,9 +25,7 @@ def fast_norm(v: Union[Tuple[float, ...], np.ndarray]) -> float:


def fast_2d_point_in_simplex(
point: Tuple[float, ...],
simplex: Union[List[Tuple[float, ...]], np.ndarray],
eps: float = 1e-8,
point: Point, simplex: SimplexPoints, eps: float = 1e-8
) -> Union[bool, np.bool_]:
(p0x, p0y), (p1x, p1y), (p2x, p2y) = simplex
px, py = point
Expand All @@ -38,9 +40,7 @@ def fast_2d_point_in_simplex(
return (t >= -eps) and (s + t <= 1 + eps)


def point_in_simplex(
point: Any, simplex: Simplex, eps: float = 1e-8
) -> Union[bool, np.bool_]:
def point_in_simplex(point: Point, simplex: SimplexPoints, eps: float = 1e-8) -> bool:
if len(point) == 2:
return fast_2d_point_in_simplex(point, simplex, eps)

Expand All @@ -51,7 +51,7 @@ def point_in_simplex(
return all(alpha > -eps) and sum(alpha) < 1 + eps


def fast_2d_circumcircle(points: np.ndarray,) -> Tuple[Tuple[float, float], float]:
def fast_2d_circumcircle(points: Iterable[Point]) -> Tuple[Tuple[float, float], float]:
"""Compute the center and radius of the circumscribed circle of a triangle
Parameters
Expand Down Expand Up @@ -88,7 +88,7 @@ def fast_2d_circumcircle(points: np.ndarray,) -> Tuple[Tuple[float, float], floa


def fast_3d_circumcircle(
points: np.ndarray,
points: Iterable[Point],
) -> Tuple[Tuple[float, float, float], float]:
"""Compute the center and radius of the circumscribed shpere of a simplex.
Expand Down Expand Up @@ -140,7 +140,7 @@ def fast_det(matrix: np.ndarray) -> float:
return np.linalg.det(matrix)


def circumsphere(pts: np.ndarray,) -> Tuple[Tuple[float, ...], float]:
def circumsphere(pts: np.ndarray) -> Tuple[Tuple[float, ...], float]:
dim = len(pts) - 1
if dim == 2:
return fast_2d_circumcircle(pts)
Expand Down Expand Up @@ -193,10 +193,12 @@ def orientation(face: np.ndarray, origin: np.ndarray) -> int:


def is_iterable_and_sized(obj: Any) -> bool:
return isinstance(obj, Iterable) and isinstance(obj, Sized)
return isinstance(obj, collections.abc.Iterable) and isinstance(
obj, collections.abc.Sized
)


def simplex_volume_in_embedding(vertices: List[Tuple[float, ...]]) -> float:
def simplex_volume_in_embedding(vertices: Iterable[Point]) -> 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.
Expand Down Expand Up @@ -277,7 +279,7 @@ class Triangulation:
or more simplices in the
"""

def __init__(self, coords: np.ndarray) -> None:
def __init__(self, coords: Iterable[Point]) -> None:
if not is_iterable_and_sized(coords):
raise TypeError("Please provide a 2-dimensional list of points")
coords = list(coords)
Expand Down Expand Up @@ -305,10 +307,10 @@ def __init__(self, coords: np.ndarray) -> None:
"(the points are linearly dependent)"
)

self.vertices = list(coords)
self.simplices = set()
self.vertices: List[Point] = list(coords)
self.simplices: Set[Simplex] = set()
# initialise empty set for each vertex
self.vertex_to_simplices = [set() for _ in coords]
self.vertex_to_simplices: List[Set[Simplex]] = [set() for _ in coords]

# find a Delaunay triangulation to start with, then we will throw it
# away and continue with our own algorithm
Expand All @@ -328,16 +330,16 @@ def add_simplex(self, simplex: Simplex) -> None:
for vertex in simplex:
self.vertex_to_simplices[vertex].add(simplex)

def get_vertices(self, indices: Sequence[int]) -> Any:
def get_vertices(self, indices: Sequence[int]) -> List[Optional[Point]]:
return [self.get_vertex(i) for i in indices]

def get_vertex(self, index: Optional[int]) -> Any:
def get_vertex(self, index: Optional[int]) -> Optional[Point]:
if index is None:
return None
return self.vertices[index]

def get_reduced_simplex(
self, point: Any, simplex: Simplex, eps: float = 1e-8
self, point: Point, simplex: Simplex, eps: float = 1e-8
) -> list:
"""Check whether vertex lies within a simplex.
Expand All @@ -364,12 +366,12 @@ def get_reduced_simplex(
return [simplex[i] for i in result]

def point_in_simplex(
self, point: Any, simplex: Simplex, eps: float = 1e-8
) -> Union[bool, np.bool_]:
self, point: Point, simplex: Simplex, eps: float = 1e-8
) -> bool:
vertices = self.get_vertices(simplex)
return point_in_simplex(point, vertices, eps)

def locate_point(self, point: Any) -> Any:
def locate_point(self, point: Point) -> Simplex:
"""Find to which simplex the point belongs.
Return indices of the simplex containing the point.
Expand All @@ -385,8 +387,11 @@ def dim(self) -> int:
return len(self.vertices[0])

def faces(
self, dim: None = None, simplices: Optional[Any] = None, vertices: None = None
) -> Iterator[Any]:
self,
dim: Optional[int] = None,
simplices: Optional[Iterable[Simplex]] = None,
vertices: Optional[Iterable[int]] = None,
) -> Iterator[Tuple[int, ...]]:
"""Iterator over faces of a simplex or vertex sequence."""
if dim is None:
dim = self.dim
Expand All @@ -407,11 +412,11 @@ def faces(
else:
return faces

def containing(self, face):
def containing(self, face: Tuple[int, ...]) -> Set[Simplex]:
"""Simplices containing a face."""
return set.intersection(*(self.vertex_to_simplices[i] for i in face))

def _extend_hull(self, new_vertex: Any, eps: float = 1e-8) -> Any:
def _extend_hull(self, new_vertex: Point, eps: float = 1e-8) -> Set[Simplex]:
# count multiplicities in order to get all hull faces
multiplicities = Counter(face for face in self.faces())
hull_faces = [face for face, count in multiplicities.items() if count == 1]
Expand Down Expand Up @@ -471,7 +476,7 @@ def circumscribed_circle(

def point_in_cicumcircle(
self, pt_index: int, simplex: Simplex, transform: np.ndarray
) -> np.bool_:
) -> bool:
# return self.fast_point_in_circumcircle(pt_index, simplex, transform)
eps = 1e-8

Expand All @@ -487,9 +492,9 @@ def default_transform(self) -> np.ndarray:
def bowyer_watson(
self,
pt_index: int,
containing_simplex: Optional[Any] = None,
containing_simplex: Optional[Simplex] = None,
transform: Optional[np.ndarray] = None,
) -> Any:
) -> Tuple[Set[Simplex], Set[Simplex]]:
"""Modified Bowyer-Watson point adding algorithm.
Create a hole in the triangulation around the new point,
Expand Down Expand Up @@ -549,7 +554,7 @@ def bowyer_watson(
new_triangles = self.vertex_to_simplices[pt_index]
return bad_triangles - new_triangles, new_triangles - bad_triangles

def _simplex_is_almost_flat(self, simplex: Simplex) -> np.bool_:
def _simplex_is_almost_flat(self, simplex: Simplex) -> bool:
return self._relative_volume(simplex) < 1e-8

def _relative_volume(self, simplex: Simplex) -> float:
Expand All @@ -565,8 +570,8 @@ def _relative_volume(self, simplex: Simplex) -> float:

def add_point(
self,
point: Any,
simplex: Optional[Any] = None,
point: Point,
simplex: Optional[Simplex] = None,
transform: Optional[np.ndarray] = None,
) -> Any:
"""Add a new vertex and create simplices as appropriate.
Expand All @@ -575,13 +580,13 @@ def add_point(
----------
point : float vector
Coordinates of the point to be added.
transform : N*N matrix of floats
Multiplication matrix to apply to the point (and neighbouring
simplices) when running the Bowyer Watson method.
simplex : tuple of ints, optional
Simplex containing the point. Empty tuple indicates points outside
the hull. If not provided, the algorithm costs O(N), so this should
be used whenever possible.
transform : N*N matrix of floats
Multiplication matrix to apply to the point (and neighbouring
simplices) when running the Bowyer Watson method.
"""
point = tuple(point)
if simplex is None:
Expand Down Expand Up @@ -626,7 +631,7 @@ def volume(self, simplex: Simplex) -> float:
def volumes(self) -> List[float]:
return [self.volume(sim) for sim in self.simplices]

def reference_invariant(self):
def reference_invariant(self) -> bool:
"""vertex_to_simplices and simplices are compatible."""
for vertex in range(len(self.vertices)):
if any(vertex not in tri for tri in self.vertex_to_simplices[vertex]):
Expand All @@ -640,26 +645,28 @@ def vertex_invariant(self, vertex):
"""Simplices originating from a vertex don't overlap."""
raise NotImplementedError

def get_neighbors_from_vertices(self, simplex: Simplex) -> Any:
def get_neighbors_from_vertices(self, simplex: Simplex) -> Set[Simplex]:
return set.union(*[self.vertex_to_simplices[p] for p in simplex])

def get_face_sharing_neighbors(self, neighbors: Any, simplex: Simplex) -> Any:
def get_face_sharing_neighbors(
self, neighbors: Set[Simplex], simplex: Simplex
) -> Set[Simplex]:
"""Keep only the simplices sharing a whole face with simplex."""
return {
simpl for simpl in neighbors if len(set(simpl) & set(simplex)) == self.dim
} # they share a face

def get_simplices_attached_to_points(self, indices: Any) -> Any:
def get_simplices_attached_to_points(self, indices: Simplex) -> Set[Simplex]:
# Get all simplices that share at least a point with the simplex
neighbors = self.get_neighbors_from_vertices(indices)
return self.get_face_sharing_neighbors(neighbors, indices)

def get_opposing_vertices(self, simplex: Simplex,) -> Any:
def get_opposing_vertices(self, simplex: Simplex) -> Tuple[int, ...]:
if simplex not in self.simplices:
raise ValueError("Provided simplex is not part of the triangulation")
neighbors = self.get_simplices_attached_to_points(simplex)

def find_opposing_vertex(vertex):
def find_opposing_vertex(vertex: int):
# find the simplex:
simp = next((x for x in neighbors if vertex not in x), None)
if simp is None:
Expand Down

0 comments on commit e959e8c

Please sign in to comment.