implement subdivide, is_simplicial_complex, is_polyhedral_fan, is_sim…

…plicial_fan; allow 3d plot

src/sage/homology/polyhedral_complex.py

@@ -27,7 +27,7 @@
 
 AUTHORS:
 
-- Yuan Zhou (2021-04): initial implementation
+- Yuan Zhou (2021-05): initial implementation
 
 List of PolyhedralComplex methods
 ---------------------------------
@@ -71,6 +71,9 @@
     :meth:`~PolyhedralComplex.is_convex` | Return ``True`` if the polyhedral complex is convex.
     :meth:`~PolyhedralComplex.is_mutable` | Return ``True`` if the polyhedral complex is mutable.
     :meth:`~PolyhedralComplex.is_immutable` | Return ``True`` if the polyhedral complex is not mutable.
+    :meth:`~PolyhedralComplex.is_simplicial_complex` | Return ``True`` if the polyhedral complex is a simplicial complex.
+    :meth:`~PolyhedralComplex.is_polyhedral_fan` | Return ``True`` if the polyhedral complex is a fan.
+    :meth:`~PolyhedralComplex.is_simplicial_fan` | Return ``True`` if the polyhedral complex is a simplicial fan.
 
 **New polyhedral complexes from old ones**
 
@@ -87,6 +90,7 @@
     :meth:`~PolyhedralComplex.product` | Return the (Cartesian) product of this polyhedral complex with another one.
     :meth:`~PolyhedralComplex.disjoint_union` | Return the disjoint union of this polyhedral complex with another one.
     :meth:`~PolyhedralComplex.join` | Return the join of this polyhedral complex with another one.
+    :meth:`~PolyhedralComplex.subdivide` | Return a new polyhedral complex (with option ``make_simplicial``) subdividing this one.
 
 **Update polyhedral complexe**
 
@@ -134,6 +138,7 @@
 from sage.rings.rational_field import QQ
 from sage.graphs.graph import Graph
 from sage.combinat.posets.posets import Poset
+from sage.misc.misc import powerset
 
 
 class PolyhedralComplex(GenericCellComplex):
@@ -672,7 +677,7 @@ def ambient_dimension(self):
 
     def plot(self, **kwds):
         """
-        Return a plot of the polyhedral complex, if it is of dim at most 2.
+        Return a plot of the polyhedral complex, if it is of dim at most 3.
 
         EXAMPLES::
 
@@ -682,7 +687,7 @@ def plot(self, **kwds):
             sage: pc.plot()
             Graphics object consisting of 10 graphics primitives
         """
-        if self.dimension() > 2:
+        if self.dimension() > 3:
             raise ValueError("Cannot plot in high dimension")
         return sum(cell.plot(**kwds) for cell in self.maximal_cell_iterator())
 
@@ -1984,6 +1989,242 @@ def remove_cell(self, cell, check=False):
         self._is_convex = None
         self._polyhedron = None
 
+    def is_simplicial_complex(self):
+        """
+        Test if this polyhedral complex is a simplicial complex.
+
+        A polyhedral complex is **simplicial** if all of its (maximal) cells
+        are simplices, i.e., every cell is a bounded polytope with `d+1`
+        vertices, where `d` is the dimension of the polytope.
+
+        EXAMPLES::
+
+            sage: p1 = Polyhedron(vertices=[(0, 0), (1, 1), (1, 2)])
+            sage: p2 = Polyhedron(rays=[(1, 0)])
+            sage: PolyhedralComplex([p1]).is_simplicial_complex()
+            True
+            sage: PolyhedralComplex([p2]).is_simplicial_complex()
+            False
+        """
+        return all(p.is_simplex() for p in self.maximal_cell_iterator())
+
+    def is_polyhedral_fan(self):
+        """
+        Test if this polyhedral complex is a polyhedral fan.
+
+        A polyhedral complex is a **fan** if all of its (maximal) cells
+        are cones.
+
+        EXAMPLES::
+
+            sage: p1 = Polyhedron(vertices=[(0, 0), (1, 1), (1, 2)])
+            sage: p2 = Polyhedron(rays=[(1, 0)])
+            sage: PolyhedralComplex([p1]).is_polyhedral_fan()
+            False
+            sage: PolyhedralComplex([p2]).is_polyhedral_fan()
+            True
+            sage: halfplane = Polyhedron(rays=[(1, 0), (-1, 0), (0, 1)])
+            sage: PolyhedralComplex([halfplane]).is_polyhedral_fan()
+            True
+        """
+        return all((p.n_vertices() == 1) and (
+                   vector(p.vertices_list()[0]) == p.ambient_space().zero())
+                   for p in self.maximal_cell_iterator())
+
+    def is_simplicial_fan(self):
+        """
+        Test if this polyhedral complex is a simplicial fan.
+
+        A polyhedral complex is a **simplicial fan** if all of its (maximal)
+        cells are simplical cones, i.e., every cell is a pointed cone (with
+        vertex being the origin) generated by `d` linearly independent rays,
+        where `d` is the dimension of the cone.
+
+        EXAMPLES::
+
+            sage: p1 = Polyhedron(vertices=[(0, 0), (1, 1), (1, 2)])
+            sage: p2 = Polyhedron(rays=[(1, 0)])
+            sage: PolyhedralComplex([p1]).is_simplicial_fan()
+            False
+            sage: PolyhedralComplex([p2]).is_simplicial_fan()
+            True
+            sage: halfplane = Polyhedron(rays=[(1, 0), (-1, 0), (0, 1)])
+            sage: PolyhedralComplex([halfplane]).is_simplicial_fan()
+            False
+        """
+        return self.is_polyhedral_fan() and all(
+               (p.n_lines() == 0 and p.n_rays() == p.dimension())
+               for p in self.maximal_cell_iterator())
+
+    def subdivide(self, make_simplicial=False,
+                  new_vertices=None, new_rays=None):
+        """
+        Construct a new polyhedral complex by iterative stellar subdivision of
+        ``self`` for each new vertex/ray given.
+
+        Currently, subdivision is only supported for bounded polyhedral complex
+        or polyhedral fan.
+
+        :param make_simplicial: boolean; optional, default ``False``.
+            If ``True``, the returned polyhedral complex is simplicial.
+
+        :param new_vertices, new_rays: list; optional, default ``None``.
+            New generators to be added during subdivision.
+
+        EXAMPLES::
+
+            sage: square_vertices = [(1, 1, 1), (-1, 1, 1), (-1, -1, 1), (1, -1, 1)]
+            sage: pc = PolyhedralComplex([
+            ....: Polyhedron(vertices=[(0, 0, 0)] + square_vertices),
+            ....: Polyhedron(vertices=[(0, 0, 2)] + square_vertices)])
+            sage: pc.is_compact() and not pc.is_simplicial_complex()
+            True
+            sage: subdivided_pc = pc.subdivide(new_vertices=[(0, 0, 1)])
+            sage: subdivided_pc
+            Polyhedral complex with 8 maximal cells
+            sage: subdivided_pc.is_simplicial_complex()
+            True
+            sage: simplicial_pc = pc.subdivide(make_simplicial=True)
+            sage: simplicial_pc
+            Polyhedral complex with 4 maximal cells
+            sage: simplicial_pc.is_simplicial_complex()
+            True
+
+            sage: fan = PolyhedralComplex([Polyhedron(rays=square_vertices)])
+            sage: fan.is_polyhedral_fan() and not fan.is_simplicial_fan()
+            True
+            sage: fan.subdivide(new_vertices=[(0, 0, 1)])
+            Traceback (most recent call last):
+            ...
+            ValueError: new vertices cannot be used for subdivision.
+            sage: subdivided_fan = fan.subdivide(new_rays=[(0, 0, 1)])
+            sage: subdivided_fan
+            Polyhedral complex with 4 maximal cells
+            sage: subdivided_fan.is_simplicial_fan()
+            True
+            sage: simplicial_fan = fan.subdivide(make_simplicial=True)
+            sage: simplicial_fan
+            Polyhedral complex with 2 maximal cells
+            sage: simplicial_fan.is_simplicial_fan()
+            True
+
+            sage: halfspace = PolyhedralComplex([Polyhedron(rays=[(0, 0, 1)],
+            ....:             lines=[(1, 0, 0), (0, 1, 0)])])
+            sage: halfspace.is_simplicial_fan()
+            False
+            sage: subdiv_halfspace = halfspace.subdivide(make_simplicial=True)
+            sage: subdiv_halfspace
+            Polyhedral complex with 4 maximal cells
+            sage: subdiv_halfspace.is_simplicial_fan()
+            True
+        """
+        if self.is_compact():
+            if new_rays:
+                raise ValueError("rays/lines cannot be used for subdivision.")
+            # bounded version of `fan.subdivide`; not require rational.
+            vertices = set([])
+            if make_simplicial and not self.is_simplicial_complex():
+                for p in self.maximal_cell_iterator():
+                    for v in p.vertices_list():
+                        vertices.add(tuple(v))
+            if new_vertices:
+                for v in new_vertices:
+                    vertices.add(tuple(v))
+            if not vertices:
+                return self    # Nothing has to be done
+            # bounded version of `fan._subdivide_stellar`; not require rational.
+            cells = list(self.maximal_cell_iterator())
+            for v in vertices:
+                new = []
+                for cell in cells:
+                    if v in cell:
+                        for cell_facet in cell.facets():
+                            facet = cell_facet.as_polyhedron()
+                            if v in facet:
+                                continue
+                            p = facet.convex_hull(Polyhedron(vertices=[v]))
+                            new.append(p)
+                    else:
+                        new.append(cell)
+                cells = new
+            return PolyhedralComplex(cells, maximality_check=False)
+        elif self.is_polyhedral_fan():
+            if new_vertices and any(vi != 0 for v in new_vertices for vi in v):
+                raise ValueError("new vertices cannot be used for subdivision.")
+            # mimic :meth:`~sage.geometry.fan <RationalPolyhedralFan>.subdivide`
+            # but here we allow for non-pointed cones, and we subdivide them.
+            rays_normalized = set([])
+            self_rays = []
+            cones = []
+            for p in self.maximal_cell_iterator():
+                prays = p.rays_list()
+                for r in prays:
+                    r_n = vector(r).normalized()
+                    r_n.set_immutable()
+                    if r_n not in rays_normalized:
+                        rays_normalized.add(r_n)
+                        self_rays.append(vector(r))
+                plines = p.lines_list()
+                if not plines:
+                    cones.append(p)
+                    continue
+                # consider a line as two rays
+                for pl in plines:
+                    l_plus = vector(pl).normalized()
+                    l_plus.set_immutable()
+                    if l_plus not in rays_normalized:
+                        rays_normalized.add(l_plus)
+                        self_rays.append(vector(pl))
+                    l_minus = (-vector(pl)).normalized()
+                    l_minus.set_immutable()
+                    if l_minus not in rays_normalized:
+                        rays_normalized.add(l_minus)
+                        self_rays.append(-vector(pl))
+                # subdivide the non-pointed p into pointed cones
+                # we rely on the canonical V-repr of Sage polyhedra.
+                num_lines = len(plines)
+                for neg_rays in powerset(range(num_lines)):
+                    lines = [vector(plines[i]) if i not in neg_rays
+                             else -vector(plines[i]) for i in range(num_lines)]
+                    cones.append(Polyhedron(rays=(prays + lines)))
+            rays = []
+            if new_rays:
+                for r in new_rays:
+                    if vector(r).is_zero():
+                        raise ValueError("zero cannot be used for subdivision.")
+                    r_n = vector(r).normalized()
+                    r_n.set_immutable()
+                    if r_n not in rays_normalized:
+                        rays_normalized.add(r_n)
+                        rays.append(vector(r))
+            if make_simplicial and not self.is_simplicial_fan():
+                rays = self_rays + rays
+            if not rays:
+                return self    # Nothing has to be done
+            # mimic :class:`RationalPolyhedralFan`._subdivide_stellar(rays)
+            # start with self maximal cells (subdivided into pointed cones)
+            for ray in rays:
+                new = []
+                for cone in cones:
+                    if ray in cone:
+                        for cone_facet in cone.facets():
+                            facet = cone_facet.as_polyhedron()
+                            if ray in facet:
+                                continue
+                            new_cone = facet.convex_hull(Polyhedron(rays=[ray]))
+                            new.append(new_cone)
+                    else:
+                        new.append(cone)
+                cones = new
+            return PolyhedralComplex(cones, maximality_check=False)
+        else:
+            # TODO: `self`` is unbounded, make it projectively simplicial.
+            # (1) homogenize self of dim d to fan in space of dim d+1;
+            # (2) call fan.subdivide(make_simplicial=True);
+            # (3) take section back to the space of dim d.
+            raise NotImplementedError('subdivision of a non-compact polyhedral ' +
+                                      'complex that is not a fan is not supported')
+
 ############################################################
 # Helper functions
 ############################################################