# Ticket #12159: trac_12159_placing_triangulation.patch

File trac_12159_placing_triangulation.patch, 39.0 KB (added by Volker Braun, 11 years ago)

Initial patch

• ## sage/geometry/polyhedra.py

# HG changeset patch
# User Volker Braun <vbraun@stp.dias.ie>
# Date 1323978224 0
# Node ID 9d0d38fe091ce0784d9ba8a4263eb1fc5836988a
Trac #12159: Placing triangulation and normal cones

This patch implements the placing triangulation.

diff --git a/sage/geometry/polyhedra.py b/sage/geometry/polyhedra.py
 a :meth:Vrepresentation objects. EXAMPLES:: sage: cube = polytopes.n_cube(3) sage: triangulation = cube.triangulate(engine='internal') # to make doctest independent of TOPCOM sage: triangulation (<0,1,2,4>, <1,2,3,4>, <1,3,4,5>, <2,3,4,6>, <3,4,5,6>, <3,5,6,7>) (<0,1,2,7>, <0,1,4,7>, <0,2,4,7>, <1,2,3,7>, <1,4,5,7>, <2,4,6,7>) sage: simplex_indices = triangulation[0]; simplex_indices (0, 1, 2, 4) (0, 1, 2, 7) sage: simplex_vertices = [ cube.Vrepresentation(i) for i in simplex_indices ] sage: simplex_vertices [A vertex at (1, 1, 1), A vertex at (-1, 1, 1), A vertex at (1, -1, 1), A vertex at (1, 1, -1)] A vertex at (1, -1, 1), A vertex at (-1, -1, -1)] sage: Polyhedron(simplex_vertices) A 3-dimensional polyhedron in QQ^3 defined as the convex hull of 4 vertices. """
• ## sage/geometry/triangulation/base.pyx

diff --git a/sage/geometry/triangulation/base.pyx b/sage/geometry/triangulation/base.pyx
 a from sage.categories.sets_cat import Sets from sage.matrix.constructor import matrix from sage.misc.misc import uniq from sage.misc.cachefunc import cached_method from functions cimport binomial from triangulations cimport \ INPUT: - point_configuration -- :class:PointConfiguration_base. The point configuration to which the point belongs. - i -- integer. The index of the point in the point configuration. EXAMPLES:: sage: pc = PointConfiguration([(0,0)]) sage: from sage.geometry.triangulation.base import Point sage: Point(123, (0,0,1), (0,0), (0,)) sage: Point(pc, 123, (0,0,1), (0,0), ()) P(0, 0) """ cdef int _index cdef tuple _projective, _affine, _reduced_affine cdef object _point_configuration cdef object _reduced_affine_vector, _reduced_projective_vector def __init__(self, i, projective, affine, reduced): def __init__(self, point_configuration, i, projective, affine, reduced): r""" Construct a :class:Point. EXAMPLES:: sage: pc = PointConfiguration([(0,0)]) sage: from sage.geometry.triangulation.base import Point sage: Point(123, (0,0,1), (0,0), (0,))   # indirect doctest sage: Point(pc, 123, (0,0,1), (0,0), ())   # indirect doctest P(0, 0) """ self._index = i self._projective = tuple(projective) self._affine = tuple(affine) self._reduced_affine = tuple(reduced) self._point_configuration = point_configuration V = point_configuration.reduced_affine_vector_space() self._reduced_affine_vector = V(self._reduced_affine) P = point_configuration.reduced_projective_vector_space() self._reduced_projective_vector = P(self.reduced_projective()) cpdef point_configuration(self): r""" Return the point configuration to which the point belongs. OUTPUT: A :class:~sage.geometry.triangulation.point_configuration.PointConfiguration. EXAMPLES: sage: pc = PointConfiguration([ (0,0), (1,0), (0,1) ]) sage: p = pc.point(0) sage: p is pc.point(0) True sage: p.point_configuration() is pc True """ return self._point_configuration def __iter__(self): EXAMPLES:: sage: pc = PointConfiguration([(0,0)]) sage: from sage.geometry.triangulation.base import Point sage: p = Point(123, (1,2,1), (3,4), (5,)) sage: p = Point(pc, 123, (1,2,1), (3,4), ()) sage: list(p)  # indirect doctest [3, 4] """ EXAMPLES:: sage: pc = PointConfiguration([(0,0)]) sage: from sage.geometry.triangulation.base import Point sage: p = Point(123, (1,2,1), (3,4), (5,)) sage: p = Point(pc, 123, (1,2,1), (3,4), ()) sage: len(p) 2 sage: p.__len__() r""" Return the projective coordinates of the point in the ambient space. EXAMPLES:: sage: pc = PointConfiguration([[10, 0, 1], [10, 0, 0], [10, 2, 3]]) sage: p = pc.point(2); p P(10, 2, 3) sage: p.affine() (10, 2, 3) sage: p.projective() (10, 2, 3, 1) sage: p.reduced_affine() (2, 2) sage: p.reduced_projective() (2, 2, 1) """ return self._projective OUTPUT: cpdef affine(self): r""" Return the affine coordinates of the point in the ambient space. A tuple containing the coordinates. EXAMPLES:: (2, 2) sage: p.reduced_projective() (2, 2, 1) sage: p.reduced_affine_vector() (2, 2) """ return self._affine return self._projective cpdef reduced_affine(self): cpdef affine(self): r""" Return the affine coordinates of the point on the hyperplane spanned by the point configuration. Return the affine coordinates of the point in the ambient space. OUTPUT: A tuple containing the coordinates. EXAMPLES:: (2, 2) sage: p.reduced_projective() (2, 2, 1) sage: p.reduced_affine_vector() (2, 2) """ return self._reduced_affine return self._affine cpdef reduced_affine(self): r""" Return the affine coordinates of the point on the hyperplane spanned by the point configuration. cpdef reduced_projective(self): r""" Return the projective coordinates of the point on the hyperplane spanned by the point configuration. OUTPUT: A tuple containing the coordinates. EXAMPLES:: (2, 2) sage: p.reduced_projective() (2, 2, 1) sage: p.reduced_affine_vector() (2, 2) """ return self._reduced_affine cpdef reduced_projective(self): r""" Return the projective coordinates of the point on the hyperplane spanned by the point configuration. OUTPUT: A tuple containing the coordinates. EXAMPLES:: sage: pc = PointConfiguration([[10, 0, 1], [10, 0, 0], [10, 2, 3]]) sage: p = pc.point(2); p P(10, 2, 3) sage: p.affine() (10, 2, 3) sage: p.projective() (10, 2, 3, 1) sage: p.reduced_affine() (2, 2) sage: p.reduced_projective() (2, 2, 1) sage: p.reduced_affine_vector() (2, 2) """ return tuple(self._reduced_affine)+(1,) cpdef reduced_affine_vector(self): """ Return the affine coordinates of the point on the hyperplane spanned by the point configuration. OUTPUT: A tuple containing the coordinates. EXAMPLES:: sage: pc = PointConfiguration([[10, 0, 1], [10, 0, 0], [10, 2, 3]]) sage: p = pc.point(2); p P(10, 2, 3) sage: p.affine() (10, 2, 3) sage: p.projective() (10, 2, 3, 1) sage: p.reduced_affine() (2, 2) sage: p.reduced_projective() (2, 2, 1) sage: p.reduced_affine_vector() (2, 2) """ return self._reduced_affine_vector cpdef reduced_projective_vector(self): """ Return the affine coordinates of the point on the hyperplane spanned by the point configuration. OUTPUT: A tuple containing the coordinates. EXAMPLES:: sage: pc = PointConfiguration([[10, 0, 1], [10, 0, 0], [10, 2, 3]]) sage: p = pc.point(2); p P(10, 2, 3) sage: p.affine() (10, 2, 3) sage: p.projective() (10, 2, 3, 1) sage: p.reduced_affine() (2, 2) sage: p.reduced_projective() (2, 2, 1) sage: p.reduced_affine_vector() (2, 2) sage: type(p.reduced_affine_vector()) """ return self._reduced_projective_vector cpdef _repr_(self): """ Return a string representation of the point. OUTPUT: String. EXAMPLES:: sage: pc = PointConfiguration([[10, 0, 1], [10, 0, 0]]) sage: from sage.geometry.triangulation.base import Point sage: p = Point(123, (0,0,1), (0,0), (0,)) sage: p = Point(pc, 123, (0,0,1), (0,0), (0,)) sage: p._repr_() 'P(0, 0)' """ :class:~sage.geometry.triangulation.point_configuration.PointConfiguration. """ def __init__(self, points): def __init__(self, points, defined_affine): r""" Construct a :class:PointConfiguration_base. - points -- a tuple of tuples of projective coordinates with 1 as the final coordinate. - defined_affine -- Boolean. Whether the point configuration is defined as a configuration of affine (as opposed to projective) points. TESTS:: sage: from sage.geometry.triangulation.base import PointConfiguration_base sage: PointConfiguration_base(((1,2,1),(2,3,1),(3,4,1)))  # indirect doctest sage: PointConfiguration_base(((1,2,1),(2,3,1),(3,4,1)), False)  # indirect doctest """ Parent.__init__(self, category = Sets()) self._init_points(points) self._is_affine = defined_affine cpdef tuple _pts cpdef int _ambient_dim cpdef int _dim cpdef object _base_ring cdef tuple _pts cdef int _ambient_dim cdef int _dim cdef object _base_ring cdef bint _is_affine cdef object _reduced_affine_vector_space, _reduced_projective_vector_space cdef _init_points(self, tuple projective_points): """ red = matrix([ red.row(i) for i in red.pivot_rows()]) else: red = matrix(0,1) self._dim = red.nrows() self._pts = tuple([ Point(i, proj.column(i), aff.column(i), red.column(i)) from sage.modules.free_module import VectorSpace self._reduced_affine_vector_space = VectorSpace(self._base_ring.fraction_field(), self._dim) self._reduced_projective_vector_space = VectorSpace(self._base_ring.fraction_field(), self._dim+1) self._pts = tuple([ Point(self, i, proj.column(i), aff.column(i), red.column(i)) for i in range(0,n) ]) self._dim = len(self._pts[0].reduced_affine()) cpdef reduced_affine_vector_space(self): """ Return the vector space that contains the affine points. OUTPUT: A vector space over the fraction field of :meth:base_ring. EXAMPLES:: sage: p = PointConfiguration([[0,0,0], [1,2,3]]) sage: p.base_ring() Integer Ring sage: p.reduced_affine_vector_space() Vector space of dimension 1 over Rational Field sage: p.reduced_projective_vector_space() Vector space of dimension 2 over Rational Field """ return self._reduced_affine_vector_space cpdef reduced_projective_vector_space(self): """ Return the vector space that is spanned by the homogeneous coordinates. OUTPUT: A vector space over the fraction field of :meth:base_ring. EXAMPLES:: sage: p = PointConfiguration([[0,0,0], [1,2,3]]) sage: p.base_ring() Integer Ring sage: p.reduced_affine_vector_space() Vector space of dimension 1 over Rational Field sage: p.reduced_projective_vector_space() Vector space of dimension 2 over Rational Field """ return self._reduced_projective_vector_space cpdef ambient_dim(self): """ Return the dimension of the ambient space of the point EXAMPLES:: sage: p = PointConfiguration([[0,0,0]]); sage: p = PointConfiguration([[0,0,0]]) sage: p.ambient_dim() 3 sage: p.dim() EXAMPLES:: sage: p = PointConfiguration([[0,0,0]]); sage: p = PointConfiguration([[0,0,0]]) sage: p.ambient_dim() 3 sage: p.dim() EXAMPLES:: sage: p = PointConfiguration([(0,0)]); sage: p = PointConfiguration([(0,0)]) sage: p.base_ring() Integer Ring sage: p = PointConfiguration([(1/2,3)]); sage: p = PointConfiguration([(1/2,3)]) sage: p.base_ring() Rational Field sage: p = PointConfiguration([(0.2, 5)]); sage: p = PointConfiguration([(0.2, 5)]) sage: p.base_ring() Real Field with 53 bits of precision """ return self._base_ring cpdef bint is_affine(self): """ Whether the configuration is defined by affine points. OUTPUT: Boolean. If true, the homogeneous coordinates all have 1 as their last entry. EXAMPLES:: sage: p = PointConfiguration([(0.2, 5), (3, 0.1)]) sage: p.is_affine() True sage: p = PointConfiguration([(0.2, 5, 1), (3, 0.1, 1)], projective=True) sage: p.is_affine() False """ return self._is_affine def _assert_is_affine(self): """ Raise a ValueError if the point configuration is not defined by affine points. EXAMPLES:: sage: p = PointConfiguration([(0.2, 5), (3, 0.1)]) sage: p._assert_is_affine() sage: p = PointConfiguration([(0.2, 5, 1), (3, 0.1, 1)], projective=True) sage: p._assert_is_affine() Traceback (most recent call last): ... ValueError: The point configuration contains projective points. """ if not self.is_affine(): raise ValueError('The point configuration contains projective points.') def __getitem__(self, i): """ Return the i-th point. EXAMPLES:: sage: p = PointConfiguration([[1,0], [2,3], [3,2]]); sage: p = PointConfiguration([[1,0], [2,3], [3,2]]) sage: [ p[i] for i in range(0,p.n_points()) ] [P(1, 0), P(2, 3), P(3, 2)] sage: list(p)
• ## sage/geometry/triangulation/element.py

diff --git a/sage/geometry/triangulation/element.py b/sage/geometry/triangulation/element.py
 a EXAMPLES:: sage: PointConfiguration.set_engine('internal')   # to make doctests independent of TOPCOM sage: p = [[0,-1,-1],[0,0,1],[0,1,0], [1,-1,-1],[1,0,1],[1,1,0]] sage: points = PointConfiguration(p) sage: triang = points.triangulate() sage: p = PointConfiguration([[0,0],[1,0],[2,1],[1,2],[0,1]]) sage: p.triangulate().gkz_phi() (3, 1, 5, 2, 4) sage: p.lexicographic_triangulation().gkz_phi() (1, 3, 4, 2, 5) NOTE: EXAMPLES:: sage: p = polytopes.cuboctahedron() sage: sc = p.triangulate().simplicial_complex() sage: sc = p.triangulate(engine='internal').simplicial_complex() sage: sc Simplicial complex with 12 vertices and 16 facets Simplicial complex with 12 vertices and 17 facets Any convex set is contractable, so its reduced homology groups vanish:: sage: sc.homology() {0: 0, 1: 0, 2: 0, 3: 0} """ from sage.homology.simplicial_complex import SimplicialComplex from sage.misc.all import flatten
• ## sage/geometry/triangulation/point_configuration.py

diff --git a/sage/geometry/triangulation/point_configuration.py b/sage/geometry/triangulation/point_configuration.py
 a Finding a single triangulation and listing all connected triangulations is implemented natively in this package. However, for more advanced options [TOPCOM]_ needs to be installed. You can find an experimental spkg at http://trac.sagemath.org/sage_trac/ticket/8169 more advanced options [TOPCOM]_ needs to be installed. It is available as an optional package for Sage, and you can install it with the command:: NOTE: sage: install_package('TOPCOM')     # not tested TOPCOM and the internal algorithms tend to enumerate triangulations in a different order. This is why we always explicitly specify the engine as engine='TOPCOM' or engine='internal' in the doctests. In your own applications, you do not need to specify the engine. By default, TOPCOM is used if it is available and the internal algorithms are used otherwise. .. note:: TOPCOM and the internal algorithms tend to enumerate triangulations in a different order. This is why we always explicitly specify the engine as engine='TOPCOM' or engine='internal' in the doctests. In your own applications, you do not need to specify the engine. By default, TOPCOM is used if it is available and the internal algorithms are used otherwise. EXAMPLES: First, we select the internal implementation for enumerating triangulations:: sage: PointConfiguration.set_engine('internal')   # to make doctests independent of TOPCOM sage: PointConfiguration.set_engine('internal')   # to make doctests independent of TOPCOM A 2-dimensional point configuration:: sage: pc.dim() 3 sage: pc.triangulate() (<0,1,2,3>, <1,2,3,4>, <2,3,4,5>, <3,4,5,6>) (<0,1,2,6>, <0,1,3,6>, <0,2,3,6>, <1,2,4,6>, <1,3,4,6>, <2,3,5,6>, <2,4,5,6>) sage: _ in pc.triangulations() True sage: len( pc.triangulations_list() ) 26 - Volker Braun: Cythonized parts of it, added a C++ implementation of the bistellar flip algorithm to enumerate all connected triangulations. - Volker Braun 2011: switched the triangulate() method to the placing triangulation (faster). """ ######################################################################## from sage.structure.unique_representation import UniqueRepresentation from sage.structure.element import Element from sage.misc.cachefunc import cached_method from sage.combinat.combination import Combinations from sage.rings.all import QQ, ZZ sage: pc1 is pc2   # indirect doctest True """ if isinstance(points, PointConfiguration): if isinstance(points, PointConfiguration_base): pc = points points = tuple( p.projective() for p in points ) projective = True defined_affine = pc.is_affine() elif projective: points = tuple( tuple(p) for p in points ) defined_affine = False else: points = tuple( tuple(p)+(1,) for p in points ) defined_affine = True if star!=None and star not in ZZ: star_point = tuple(star) if len(star_point)d_max: p_max = p return p_max def contained_simplex(self, large=True, initial_point=None): """ Return a simplex contained in the point configuration. INPUT: - large -- boolean. Whether to attempt to return a large simplex. - initial_point -- a :class:~sage.geometry.triangulation.base.Point or None (default). A specific point to start with when picking the simplex vertices. OUTPUT: A tuple of points that span a simplex of dimension :meth:dim. If large==True, the simplex is constructed by sucessively picking the farthest point. This will ensure that the simplex is not unneccessarily small, but will in general not return a maximal simplex. EXAMPLES:: sage: pc = PointConfiguration([(0,0),(1,0),(2,1),(1,1),(0,1)]) sage: pc.contained_simplex() (P(0, 1), P(2, 1), P(1, 0)) sage: pc.contained_simplex(large=False) (P(0, 1), P(1, 1), P(1, 0)) sage: pc.contained_simplex(initial_point=pc.point(0)) (P(0, 0), P(1, 1), P(1, 0)) sage: pc = PointConfiguration([[0,0],[0,1],[1,0],[1,1],[-1,-1]]) sage: pc.contained_simplex() (P(-1, -1), P(1, 1), P(0, 1)) TESTS:: sage: pc = PointConfiguration([[0,0],[0,1],[1,0]]) sage: pc.contained_simplex() (P(1, 0), P(0, 1), P(0, 0)) sage: pc = PointConfiguration([[0,0],[0,1]]) sage: pc.contained_simplex() (P(0, 1), P(0, 0)) sage: pc = PointConfiguration([[0,0]]) sage: pc.contained_simplex() (P(0, 0),) sage: pc = PointConfiguration([]) sage: pc.contained_simplex() () """ self._assert_is_affine() if self.n_points()==0: return tuple() points = list(self.points()) if initial_point is None: origin = points.pop() else: origin = initial_point points.remove(origin) vertices = [origin] edges = [] while len(vertices) <= self.dim(): if large: p = self.farthest_point(vertices, points) points.remove(p) else: p = points.pop() edge = p.reduced_affine_vector()-origin.reduced_affine_vector() if len(edges)>0 and (ker * edge).is_zero(): continue vertices.append(p) edges.append(edge) ker = matrix(edges).right_kernel().matrix() return tuple(vertices) def placing_triangulation(self, point_order=None): r""" Construct the placing (pushing) triangulation. INPUT: - point_order -- list of points or integers. The order in which the points are to be placed. OUTPUT: A :class:~sage.geometry.triangulation.triangulation.Triangulation. EXAMPLES:: sage: pc = PointConfiguration([(0,0),(1,0),(2,1),(1,2),(0,1)]) sage: pc.placing_triangulation() (<0,1,2>, <0,2,4>, <2,3,4>) sage: U=matrix([ ...      [ 0, 0, 0, 0, 0, 2, 4,-1, 1, 1, 0, 0, 1, 0], ...      [ 0, 0, 0, 1, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0], ...      [ 0, 2, 0, 0, 0, 0,-1, 0, 1, 0, 1, 0, 0, 1], ...      [ 0, 1, 1, 0, 0, 1, 0,-2, 1, 0, 0,-1, 1, 1], ...      [ 0, 0, 0, 0, 1, 0,-1, 0, 0, 0, 0, 0, 0, 0] ...   ]) sage: p = PointConfiguration(U.columns()) sage: triangulation = p.placing_triangulation();  triangulation (<0,2,3,4,6,7>, <0,2,3,4,6,12>, <0,2,3,4,7,13>, <0,2,3,4,12,13>, <0,2,3,6,7,13>, <0,2,3,6,12,13>, <0,2,4,6,7,13>, <0,2,4,6,12,13>, <0,3,4,6,7,12>, <0,3,4,7,12,13>, <0,3,6,7,12,13>, <0,4,6,7,12,13>, <1,3,4,5,6,12>, <1,3,4,6,11,12>, <1,3,4,7,11,13>, <1,3,4,11,12,13>, <1,3,6,7,11,13>, <1,3,6,11,12,13>, <1,4,6,7,11,13>, <1,4,6,11,12,13>, <3,4,6,7,11,12>, <3,4,7,11,12,13>, <3,6,7,11,12,13>, <4,6,7,11,12,13>) sage: sum(p.volume(t) for t in triangulation) 42 """ facet_normals = dict() def facets_of_simplex(simplex): """ Return the facets of the simplex and store the normals in facet_normals """ simplex = list(simplex) origin = simplex[0] rest = simplex[1:] span = matrix([ origin.reduced_affine_vector()-p.reduced_affine_vector() for p in rest ]) # span.inverse() linearly transforms the simplex into the unit simplex normals = span.inverse().columns() facets = [] # The facets incident to the chosen vertex "origin" for opposing_vertex, normal in zip(rest, normals): facet = frozenset([origin] + [ p for p in rest if p is not opposing_vertex ]) facets.append(facet) normal.set_immutable() facet_normals[facet] = normal # The remaining facet that is not incident to "origin" facet = frozenset(rest) normal = -sum(normals) normal.set_immutable() facet_normals[facet] = normal facets.append(facet) return set(facets) # input verification self._assert_is_affine() if point_order is None: point_order = list(self.points()) elif isinstance(point_order[0], Point): point_order = list(point_order) assert all(p.point_configuration()==self for p in point_order) else: point_order = [ self.point(i) for i in point_order ] assert all(p in self.points() for p in point_order) # construct the initial simplex simplices = [ frozenset(self.contained_simplex()) ] for s in simplices[0]: try: point_order.remove(s) except ValueError: pass facets = facets_of_simplex(simplices[0]) # successively place the remaining points for point in point_order: # identify visible facets visible_facets = [] for facet in facets: origin = iter(facet).next() normal = facet_normals[facet] v = point.reduced_affine_vector() - origin.reduced_affine_vector() if v*normal>0: visible_facets.append(facet) # construct simplices over each visible facet new_facets = set() for facet in visible_facets: simplex = frozenset(list(facet) + [point]) # print 'simplex', simplex simplices.append(simplex) for facet in facets_of_simplex(simplex): if facet in visible_facets: continue if facet in new_facets: new_facets.remove(facet) continue new_facets.add(facet) facets.difference_update(visible_facets) facets.update(new_facets) # construct the triangulation triangulation = [ [p.index() for p in simplex] for simplex in simplices ] return self(triangulation) pushing_triangulation = placing_triangulation