from collections import Counter

 from collections.abc import Iterable, Sized
 from itertools import chain, combinations
 from math import factorial

 import scipy.spatial

 from numpy import square, zeros, subtract, array, ones, dot, asarray, concatenate, average, eye, mean, abs, sqrt
 from numpy import sum as nsum
 from numpy.linalg import det as ndet
 from numpy.linalg import slogdet, solve, matrix_rank, norm
 

 def fast_norm(v):

     """ Manually take the vector norm for len 2, 3 vectors. Defaults to a square root of the dot product
     for larger vectors. 
     
     Note that for large vectors, it is possible for integer overflow to occur. 
     For instance:
     vec = [49024, 59454, 12599, -63721, 18517, 27961]
     dot(vec, vec) = -1602973744
     
     """
     len_v = len(v)

     # notice this method can be even more optimised

     if len_v == 2:

         return sqrt(v[0] * v[0] + v[1] * v[1])

     if len_v == 3:

         return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2])
     return sqrt(dot(v, v))

 def fast_2d_point_in_simplex(point, simplex, eps=1e-8):

     (p0x, p0y), (p1x, p1y), (p2x, p2y) = simplex

     px, py = point
 

     area = 0.5 * (-p1y * p2x + p0y * (p2x - p1x) + p1x * p2y + p0x * (p1y - p2y))

     s = 1 / (2 * area) * (+p0y * p2x + (p2y - p0y) * px - p0x * p2y + (p0x - p2x) * py)

     if s < -eps or s > 1 + eps:

         return False

     t = 1 / (2 * area) * (+p0x * p1y + (p0y - p1y) * px - p0y * p1x + (p1x - p0x) * py)

     return (t >= -eps) and (s + t <= 1 + eps)

 def point_in_simplex(point, simplex, eps=1e-8):

     # simplex is list

     if len(point) == 2:
         return fast_2d_point_in_simplex(point, simplex, eps)
 

     x0 = array(simplex[0], dtype=float)
     vectors = array(simplex[1:], dtype=float) - x0
     alpha = solve(vectors.T, point - x0)

 
     return all(alpha > -eps) and sum(alpha) < 1 + eps
 
 

 def fast_2d_circumcircle(points):

     """Compute the center and radius of the circumscribed circle of a triangle

 
     Parameters
     ----------

     points: 2D array-like
         the points of the triangle to investigate

 
     Returns
     -------

     tuple

         (center point : tuple(int), radius: int)

     """

     points = array(points)

     # transform to relative coordinates

     pts = points[1:] - points[0]

 
     (x1, y1), (x2, y2) = pts

     # compute the length squared

     l1 = x1 * x1 + y1 * y1
     l2 = x2 * x2 + y2 * y2

     # compute some determinants

     dx = +l1 * y2 - l2 * y1
     dy = -l1 * x2 + l2 * x1
     aa = +x1 * y2 - x2 * y1

     a = 2 * aa
 

     # compute center

     x = dx / a
     y = dy / a

     radius = sqrt(x * x + y * y)  # radius = norm([x, y])

     return (x + points[0][0], y + points[0][1]), radius

 
 
 def fast_3d_circumcircle(points):

     """Compute the center and radius of the circumscribed shpere of a simplex.

 
     Parameters
     ----------

     points: 2D array-like
         the points of the triangle to investigate

 
     Returns
     -------

     tuple

         (center point : tuple(int), radius: int)

     """

     points = array(points)

     pts = points[1:] - points[0]

 
     (x1, y1, z1), (x2, y2, z2), (x3, y3, z3) = pts

     l1 = x1 * x1 + y1 * y1 + z1 * z1
     l2 = x2 * x2 + y2 * y2 + z2 * z2
     l3 = x3 * x3 + y3 * y3 + z3 * z3
 

     # Compute some determinants:

     dx = +l1 * (y2 * z3 - z2 * y3) - l2 * (y1 * z3 - z1 * y3) + l3 * (y1 * z2 - z1 * y2)
     dy = +l1 * (x2 * z3 - z2 * x3) - l2 * (x1 * z3 - z1 * x3) + l3 * (x1 * z2 - z1 * x2)
     dz = +l1 * (x2 * y3 - y2 * x3) - l2 * (x1 * y3 - y1 * x3) + l3 * (x1 * y2 - y1 * x2)
     aa = +x1 * (y2 * z3 - z2 * y3) - x2 * (y1 * z3 - z1 * y3) + x3 * (y1 * z2 - z1 * y2)

     a = 2 * aa

     center = (dx / a, -dy / a, dz / a)

     radius = fast_norm(center)

     center = (
         center[0] + points[0][0],
         center[1] + points[0][1],
         center[2] + points[0][2],
     )

     return center, radius

 def fast_det(matrix):

     matrix = asarray(matrix, dtype=float)

     if matrix.shape == (2, 2):
         return matrix[0][0] * matrix[1][1] - matrix[1][0] * matrix[0][1]
     elif matrix.shape == (3, 3):
         a, b, c, d, e, f, g, h, i = matrix.ravel()

         return a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g)

     else:

         return ndet(matrix)

 def circumsphere(pts):
     dim = len(pts) - 1
     if dim == 2:
         return fast_2d_circumcircle(pts)
     if dim == 3:
         return fast_3d_circumcircle(pts)
 
     # Modified method from http://mathworld.wolfram.com/Circumsphere.html

     mat = array([[nsum(square(pt)), *pt, 1] for pt in pts])
     center = zeros(dim)
     a = 1 / (2 * ndet(mat[:, 1:]))
     factor = a
     ind = ones((dim + 2,), bool)

     for i in range(1, len(pts)):

         ind[i - 1] = True
         ind[i] = False
         center[i - 1] = factor * ndet(mat[:, ind])
         factor *= -1

 
     x0 = pts[0]

     vec = subtract(center, x0)
     radius = sqrt(dot(vec, vec))

 
     return tuple(center), radius
 
 

 def orientation(face, origin):

     """Compute the orientation of the face with respect to a point, origin.

 
     Parameters
     ----------

     face : array-like, of shape (N-dim, N-dim)

         The hyperplane we want to know the orientation of
         Do notice that the order in which you provide the points is critical

     origin : array-like, point of shape (N-dim)

         The point to compute the orientation from
 
     Returns
     -------

     0 if the origin lies in the same hyperplane as face,
     -1 or 1 to indicate left or right orientation

     If two points lie on the same side of the face, the orientation will
     be equal, if they lie on the other side of the face, it will be negated.

     """

     vectors = array(face)
     sign, logdet = slogdet(vectors - origin)

     if logdet < -50:  # assume it to be zero when it's close to zero
         return 0
     return sign

 def is_iterable_and_sized(obj):
     return isinstance(obj, Iterable) and isinstance(obj, Sized)

 def simplex_volume_in_embedding(vertices) -> 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.
 

     This algorithm has not been tested for numerical stability.
 

     Parameters
     ----------

     vertices : 2D arraylike of floats

 
     Returns
     -------

     volume : int

         the volume of the simplex with given vertices.
 
     Raises
     ------

     ValueError

         if the vertices do not form a simplex (for example,
         because they are coplanar, colinear or coincident).
     """
     # Implements http://mathworld.wolfram.com/Cayley-MengerDeterminant.html
     # Modified from https://codereview.stackexchange.com/questions/77593/calculating-the-volume-of-a-tetrahedron
 

     vertices = asarray(vertices, dtype=float)

     dim = len(vertices[0])
     if dim == 2:
         # Heron's formula

         a, b, c = scipy.spatial.distance.pdist(vertices, metric="euclidean")

         s = 0.5 * (a + b + c)

         return sqrt(s * (s - a) * (s - b) * (s - c))

 
     # β_ij = |v_i - v_k|²

     sq_dists = scipy.spatial.distance.pdist(vertices, metric="sqeuclidean")

 
     # Add border while compressed
     num_verts = scipy.spatial.distance.num_obs_y(sq_dists)

     bordered = concatenate((ones(num_verts), sq_dists))

 
     # Make matrix and find volume
     sq_dists_mat = scipy.spatial.distance.squareform(bordered)
 

     coeff = -((-2) ** (num_verts - 1)) * factorial(num_verts - 1) ** 2

     vol_square = fast_det(sq_dists_mat) / coeff

     if vol_square < 0:
         if abs(vol_square) < 1e-15:
             return 0

         raise ValueError("Provided vertices do not form a simplex")

     return sqrt(vol_square)

 class Triangulation:

     """A triangulation object.

     Parameters
     ----------
     coords : 2d array-like of floats

         Coordinates of vertices.

     Attributes
     ----------
     vertices : list of float tuples
         Coordinates of the triangulation vertices.
     simplices : set of integer tuples
         List with indices of vertices forming individual simplices

     vertex_to_simplices : list of sets
         Set of simplices connected to a vertex, the index of the vertex is the
         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 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 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):
             raise ValueError("Coordinates dimension mismatch")
 

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

         if len(coords) < dim + 1:

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

         coords = list(map(tuple, coords))

         vectors = subtract(coords[1:], coords[0])
         if 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)

 
     def delete_simplex(self, simplex):
         simplex = tuple(sorted(simplex))
         self.simplices.remove(simplex)
         for vertex in simplex:
             self.vertex_to_simplices[vertex].remove(simplex)
 
     def add_simplex(self, simplex):
         simplex = tuple(sorted(simplex))
         self.simplices.add(simplex)
         for vertex in simplex:
             self.vertex_to_simplices[vertex].add(simplex)
 

     def get_vertices(self, indices):

         return [self.get_vertex(i) for i in indices]
 
     def get_vertex(self, index):
         if index is None:
             return None
         return self.vertices[index]

     def get_reduced_simplex(self, point, simplex, eps=1e-8) -> list:

         """Check whether vertex lies within a simplex.
 
         Returns
         -------
         vertices : list of ints

             Indices of vertices of the simplex to which the vertex belongs.
             An empty list indicates that the vertex is outside the simplex.

         """

         # XXX: in the end we want to lose this method

         if len(simplex) != self.dim + 1:
             # We are checking whether point belongs to a face.
             simplex = self.containing(simplex).pop()

         x0 = array(self.vertices[simplex[0]])
         vectors = array(self.get_vertices(simplex[1:])) - x0
         alpha = solve(vectors.T, point - x0)

         if any(alpha < -eps) or sum(alpha) > 1 + eps:
             return []

 
         result = [i for i, a in enumerate(alpha, 1) if a > eps]
         if sum(alpha) < 1 - eps:
             result.insert(0, 0)
 

         return [simplex[i] for i in result]
 

     def point_in_simplex(self, point, simplex, eps=1e-8):

         vertices = self.get_vertices(simplex)
         return point_in_simplex(point, vertices, eps)

     def locate_point(self, point):
         """Find to which simplex the point belongs.
 
         Return indices of the simplex containing the point.
         Empty tuple means the point is outside the triangulation
         """
         for simplex in self.simplices:

             if self.point_in_simplex(point, simplex):

                 return simplex

         return ()
 
     @property
     def dim(self):
         return len(self.vertices[0])
 
     def faces(self, dim=None, simplices=None, vertices=None):
         """Iterator over faces of a simplex or vertex sequence."""
         if dim is None:
             dim = self.dim
 
         if simplices is not None and vertices is not None:
             raise ValueError("Only one of simplices and vertices is allowed.")
         if vertices is not None:
             vertices = set(vertices)
             simplices = chain(*(self.vertex_to_simplices[i] for i in vertices))
             simplices = set(simplices)
         elif simplices is None:
             simplices = self.simplices
 

         faces = (face for tri in simplices for face in combinations(tri, dim))

 
         if vertices is not None:
             return (face for face in faces if all(i in vertices for i in face))
         else:
             return faces
 
     def containing(self, face):
         """Simplices containing a face."""
         return set.intersection(*(self.vertex_to_simplices[i] for i in face))
 

     def _extend_hull(self, new_vertex, eps=1e-8):

         # 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]

         # compute the center of the convex hull, this center lies in the hull
         # we do not really need the center, we only need a point that is
         # guaranteed to lie strictly within the hull

         hull_points = self.get_vertices(self.hull)

         pt_center = average(hull_points, axis=0)

         pt_index = len(self.vertices)
         self.vertices.append(new_vertex)

 
         new_simplices = set()

         for face in hull_faces:

             # do orientation check, if orientation is the same, it lies on
             # the same side of the face, otherwise, it lies on the other
             # side of the face

             pts_face = tuple(self.get_vertices(face))

             orientation_inside = orientation(pts_face, pt_center)
             orientation_new_point = orientation(pts_face, new_vertex)
             if orientation_inside == -orientation_new_point:
                 # if the orientation of the new vertex is zero or directed
                 # towards the center, do not add the simplex

                 simplex = (*face, pt_index)
                 if not self._simplex_is_almost_flat(simplex):
                     self.add_simplex(simplex)
                     new_simplices.add(simplex)

         if len(new_simplices) == 0:

             # We tried to add an internal point, revert and raise.
             for tri in self.vertex_to_simplices[pt_index]:
                 self.simplices.remove(tri)
             del self.vertex_to_simplices[pt_index]
             del self.vertices[pt_index]
             raise ValueError("Candidate vertex is inside the hull.")
 

         return new_simplices

     def circumscribed_circle(self, simplex, transform):

         """Compute the center and radius of the circumscribed circle of a simplex.

 
         Parameters
         ----------
         simplex : tuple of ints
             the simplex to investigate
 
         Returns
         -------

         tuple (center point, radius)

             The center and radius of the circumscribed circle

         """

         pts = dot(self.get_vertices(simplex), transform)

         return circumsphere(pts)

     def point_in_cicumcircle(self, pt_index, simplex, transform):
         # return self.fast_point_in_circumcircle(pt_index, simplex, transform)

         eps = 1e-8
 

         center, radius = self.circumscribed_circle(simplex, transform)

         pt = dot(self.get_vertices([pt_index]), transform)[0]

         return norm(center - pt) < (radius * (1 + eps))

 
     @property

     def default_transform(self):

         return eye(self.dim)

     def bowyer_watson(self, pt_index, containing_simplex=None, transform=None):

         """Modified Bowyer-Watson point adding algorithm.

         Create a hole in the triangulation around the new point,
         then retriangulate this hole.

         Parameters
         ----------

         pt_index: number
             the index of the point to inspect

 
         Returns
         -------

         deleted_simplices : set of tuples
             Simplices that have been deleted
         new_simplices : set of tuples
             Simplices that have been added

         """

         queue = set()
         done_simplices = set()

         transform = self.default_transform if transform is None else transform

         if containing_simplex is None:
             queue.update(self.vertex_to_simplices[pt_index])
         else:
             queue.add(containing_simplex)
 

         bad_triangles = set()
 
         while len(queue):
             simplex = queue.pop()
             done_simplices.add(simplex)
 

             if self.point_in_cicumcircle(pt_index, simplex, transform):

                 self.delete_simplex(simplex)

                 todo_points = set(simplex)

                 bad_triangles.add(simplex)
 

                 # Get all simplices that share at least a point with the simplex

                 neighbors = self.get_neighbors_from_vertices(todo_points)

                 # Filter out the already evaluated simplices

                 neighbors = neighbors - done_simplices
                 neighbors = self.get_face_sharing_neighbors(neighbors, simplex)
                 queue.update(neighbors)

 
         faces = list(self.faces(simplices=bad_triangles))
 
         multiplicities = Counter(face for face in faces)

         hole_faces = [face for face in faces if multiplicities[face] < 2]

 
         for face in hole_faces:
             if pt_index not in face:

                 simplex = (*face, pt_index)
                 if not self._simplex_is_almost_flat(simplex):
                     self.add_simplex(simplex)

         new_triangles = self.vertex_to_simplices[pt_index]
         return bad_triangles - new_triangles, new_triangles - bad_triangles

     def _simplex_is_almost_flat(self, simplex):
         return self._relative_volume(simplex) < 1e-8
 
     def _relative_volume(self, simplex):
         """Compute the volume of a simplex divided by the average (Manhattan)
         distance of its vertices. The advantage of this is that the relative
         volume is only dependent on the shape of the simplex and not on the
         absolute size. Due to the weird scaling, the only use of this method
         is to check that a simplex is almost flat."""

         vertices = array(self.get_vertices(simplex))

         vectors = vertices[1:] - vertices[0]

         average_edge_length = mean(abs(vectors))

         return self.volume(simplex) / (average_edge_length ** self.dim)
 

     def add_point(self, point, simplex=None, transform=None):

         """Add a new vertex and create simplices as appropriate.
 
         Parameters
         ----------
         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.
         """

         point = tuple(point)

         if simplex is None:
             simplex = self.locate_point(point)
 

         actual_simplex = simplex

         self.vertex_to_simplices.append(set())

         if not simplex:

             temporary_simplices = self._extend_hull(point)

 
             pt_index = len(self.vertices) - 1

             deleted_simplices, added_simplices = self.bowyer_watson(
                 pt_index, transform=transform
             )

 
             deleted = deleted_simplices - temporary_simplices
             added = added_simplices | (temporary_simplices - deleted_simplices)
             return deleted, added

         else:

             reduced_simplex = self.get_reduced_simplex(point, simplex)

             if not reduced_simplex:

                 self.vertex_to_simplices.pop()  # revert adding vertex

                 raise ValueError("Point lies outside of the specified simplex.")

             else:
                 simplex = reduced_simplex
 
         if len(simplex) == 1:

             self.vertex_to_simplices.pop()  # revert adding vertex

             raise ValueError("Point already in triangulation.")
         else:

             pt_index = len(self.vertices)
             self.vertices.append(point)

             return self.bowyer_watson(pt_index, actual_simplex, transform)

 
     def volume(self, simplex):

         prefactor = factorial(self.dim)
         vertices = array(self.get_vertices(simplex))

         vectors = vertices[1:] - vertices[0]

         return float(abs(fast_det(vectors)) / prefactor)

     def volumes(self):
         return [self.volume(sim) for sim in self.simplices]
 

     def reference_invariant(self):
         """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]):

                 return False
         for simplex in self.simplices:

             if any(simplex not in self.vertex_to_simplices[pt] for pt in simplex):

                 return False
         return True
 
     def vertex_invariant(self, vertex):
         """Simplices originating from a vertex don't overlap."""
         raise NotImplementedError
 

     def get_neighbors_from_vertices(self, simplex):

         return set.union(*[self.vertex_to_simplices[p] for p in simplex])

 
     def get_face_sharing_neighbors(self, neighbors, 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):
         # 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):
         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):
             # find the simplex:
             simp = next((x for x in neighbors if vertex not in x), None)
             if simp is None:
                 return None
             opposing = set(simp) - set(simplex)
             assert len(opposing) == 1
             return opposing.pop()

         result = tuple(find_opposing_vertex(v) for v in simplex)
         return result
 

     @property
     def hull(self):
         """Compute hull from triangulation.

 
         Parameters
         ----------

         check : bool, default: True

             Whether to raise an error if the computed hull is different from
             stored.
 
         Returns
         -------
         hull : set of int
             Vertices in the hull.
         """
         counts = Counter(self.faces())
         if any(i > 2 for i in counts.values()):

             raise RuntimeError(
                 "Broken triangulation, a (N-1)-dimensional"
                 " appears in more than 2 simplices."
             )

         hull = {point for face, count in counts.items() if count == 1 for point in face}

         return hull
 
     def convex_invariant(self, vertex):
         """Hull is convex."""
         raise NotImplementedError