Trac #34326: ConvexSet_base.representative_point, Polyhedron_base.an_…

…affine_basis for unbounded polyhedra

`ConvexSet_base.representative_point` generalizes
`Polyhedron_base.representative_point`.
We implement it via `an_affine_basis`.

We extend the implementation of that so it handles unbounded polyhedra.

URL: https://trac.sagemath.org/34326
Reported by: mkoeppe
Ticket author(s): Matthias Koeppe
Reviewer(s): Travis Scrimshaw

src/sage/geometry/cone.py

@@ -3533,21 +3533,15 @@ def an_affine_basis(self):
 
             sage: quadrant = Cone([(1,0), (0,1)])
             sage: quadrant.an_affine_basis()
-            Traceback (most recent call last):
-            ...
-            NotImplementedError: this function is not implemented for unbounded polyhedra
+            [(0, 0), (1, 0), (0, 1)]
             sage: ray = Cone([(1, 1)])
             sage: ray.an_affine_basis()
-            Traceback (most recent call last):
-            ...
-            NotImplementedError: this function is not implemented for unbounded polyhedra
+            [(0, 0), (1, 1)]
             sage: line = Cone([(1,0), (-1,0)])
             sage: line.an_affine_basis()
-            Traceback (most recent call last):
-            ...
-            NotImplementedError: this function is not implemented for unbounded polyhedra
+            [(1, 0), (0, 0)]
         """
-        return self.polyhedron().an_affine_basis()
+        return [vector(v) for v in self.polyhedron().an_affine_basis()]
 
     @cached_method
     def strict_quotient(self):

src/sage/geometry/convex_set.py

@@ -622,6 +622,24 @@ def relative_interior(self):
             return self
         raise NotImplementedError
 
+    @cached_method
+    def representative_point(self):
+        """
+        Return a "generic" point of ``self``.
+
+        OUTPUT:
+
+        A point in the relative interior of ``self`` as a coordinate vector.
+
+        EXAMPLES::
+
+            sage: C = Cone([[1, 2, 0], [2, 1, 0]])
+            sage: C.representative_point()
+            (1, 1, 0)
+        """
+        affine_basis = self.an_affine_basis()
+        return sum(affine_basis) / len(affine_basis)
+
     def _test_convex_set(self, tester=None, **options):
         """
         Run some tests on the methods of :class:`ConvexSet_base`.
@@ -746,7 +764,6 @@ def some_elements(self):
             raise NotImplementedError
         return list(self._some_elements_())
 
-    @abstract_method(optional=True)
     def _some_elements_(self):
         r"""
         Generate some points of ``self``.
@@ -757,11 +774,12 @@ def _some_elements_(self):
 
             sage: from sage.geometry.convex_set import ConvexSet_base
             sage: C = ConvexSet_base()
-            sage: C._some_elements_(C)
+            sage: list(C._some_elements_())
             Traceback (most recent call last):
             ...
             TypeError: 'NotImplementedType' object is not callable
         """
+        yield self.representative_point()
 
     @abstract_method(optional=True)
     def cartesian_product(self, other):

src/sage/geometry/polyhedron/base1.py

@@ -426,7 +426,7 @@ def ambient_dim(self):
 
     def an_affine_basis(self):
         """
-        Return points in ``self`` that are a basis for the affine span of the polytope.
+        Return points in ``self`` that form a basis for the affine span of ``self``.
 
         This implementation of the method :meth:`~sage.geometry.convex_set.ConvexSet_base.an_affine_basis`
         for polytopes guarantees the following:
@@ -438,6 +438,9 @@ def an_affine_basis(self):
           one vertex that is in the difference. Thus this method
           is independent of the concrete realization of the polytope.
 
+        For unbounded polyhedra, the result may contain arbitrary points of ``self``,
+        not just vertices.
+
         EXAMPLES::
 
             sage: P = polytopes.cube()
@@ -455,42 +458,50 @@ def an_affine_basis(self):
              A vertex at (4, 1, 3, 5, 2),
              A vertex at (4, 2, 5, 3, 1)]
 
-        The method is not implemented for unbounded polyhedra::
+        Unbounded polyhedra::
 
-            sage: p = Polyhedron(vertices=[(0,0)],rays=[(1,0),(0,1)])
+            sage: p = Polyhedron(vertices=[(0, 0)], rays=[(1,0), (0,1)])
+            sage: p.an_affine_basis()
+            [A vertex at (0, 0), (1, 0), (0, 1)]
+            sage: p = Polyhedron(vertices=[(2, 1)], rays=[(1,0), (0,1)])
             sage: p.an_affine_basis()
-            Traceback (most recent call last):
-            ...
-            NotImplementedError: this function is not implemented for unbounded polyhedra
+            [A vertex at (2, 1), (3, 1), (2, 2)]
+            sage: p = Polyhedron(vertices=[(2, 1)], rays=[(1,0)], lines=[(0,1)])
+            sage: p.an_affine_basis()
+            [(2, 1), A vertex at (2, 0), (3, 0)]
         """
-        if not self.is_compact():
-            raise NotImplementedError("this function is not implemented for unbounded polyhedra")
-
         chain = self.a_maximal_chain()[1:]  # we exclude the empty face
         chain_indices = [face.ambient_V_indices() for face in chain]
         basis_indices = []
 
         # We use in the following that elements in ``chain_indices`` are sorted lists
         # of V-indices.
-        # Thus for each two faces we can easily find the first vertex that differs.
+        # Thus for each two faces we can easily find the first Vrep that differs.
         for dim, face in enumerate(chain_indices):
             if dim == 0:
-                # Append the vertex.
-                basis_indices.append(face[0])
+                # Append the V-indices of the minimal face.
+                basis_indices.extend(face[:])
                 continue
 
             prev_face = chain_indices[dim-1]
             for i in range(len(prev_face)):
                 if prev_face[i] != face[i]:
-                    # We found a vertex that ``face`` has, but its facet does not.
+                    # We found a Vrep that ``face`` has, but its facet does not.
                     basis_indices.append(face[i])
                     break
             else:  # no break
                 # ``prev_face`` contains all the same vertices as ``face`` until now.
-                # But ``face`` is guaranteed to contain one more vertex (at least).
+                # But ``face`` is guaranteed to contain one more Vrep (at least).
                 basis_indices.append(face[len(prev_face)])
 
-        return [self.Vrepresentation()[i] for i in basis_indices]
+        Vreps = [self.Vrepresentation()[i] for i in basis_indices]
+        if all(vrep.is_vertex() for vrep in Vreps):
+            return Vreps
+        for vrep in Vreps:
+            if vrep.is_vertex():
+                vertex = vrep
+        return [vrep if vrep.is_vertex() else vertex.vector() + vrep.vector()
+                for vrep in Vreps]
 
     @abstract_method
     def a_maximal_chain(self):

src/sage/geometry/relative_interior.py

@@ -276,6 +276,22 @@ def _some_elements_(self):
             if p in self:
                 yield p
 
+    def representative_point(self):
+        """
+        Return a "generic" point of ``self``.
+
+        OUTPUT:
+
+        A point in ``self`` (thus, in the relative interior of ``self``) as a coordinate vector.
+
+        EXAMPLES::
+
+            sage: C = Cone([[1, 2, 0], [2, 1, 0]])
+            sage: C.relative_interior().representative_point()
+            (1, 1, 0)
+        """
+        return self._polyhedron.representative_point()
+
     def _repr_(self):
         r"""
         Return a description of ``self``.