Trac #28766: Implement incidence matrix for cones

This ticket implements the incidence matrix for cones. This method is
used to simplify getting the face lattice.

Columns correspond to facets, rows to rays. So there is no extra row for
the vertex.

As a follow up, one can then obtain the combinatorial polyhedron of a
cone.

URL: https://trac.sagemath.org/28766
Reported by: gh-kliem
Ticket author(s): Jonathan Kliem
Reviewer(s): Laith Rastanawi

src/sage/geometry/cone.py

@@ -2354,9 +2354,6 @@ def face_lattice(self):
                 S = self.linear_subspace()
                 subspace_rays = []
                 atom_to_ray = []
-                atom_to_facets = []
-                normals = self.facet_normals()
-                facet_to_atoms = [[] for normal in normals]
                 for i, ray in enumerate(self):
                     # This try...except tests whether ray lies in S;
                     # "ray in S" does not work because ray lies in a
@@ -2367,12 +2364,6 @@ def face_lattice(self):
                         S(ray)
                         subspace_rays.append(i)
                     except (TypeError, ValueError):
-                        facets = [j for j, normal in enumerate(normals)
-                                    if ray * normal == 0]
-                        atom_to_facets.append(facets)
-                        atom = len(atom_to_ray)
-                        for j in facets:
-                            facet_to_atoms[j].append(atom)
                         atom_to_ray.append(i)
 
                 def ConeFace(atoms, facets):
@@ -2389,6 +2380,15 @@ def ConeFace(atoms, facets):
                     else:
                         return self
 
+                # Obtain a modified version of the incidence matrix,
+                # with rows corresponding to rays in subspace removed.
+                mod_incidence_matrix = self.incidence_matrix()[atom_to_ray]
+
+                atom_to_facets = [row.nonzero_positions()
+                                  for row in mod_incidence_matrix.rows()]
+                facet_to_atoms = [column.nonzero_positions()
+                                  for column in mod_incidence_matrix.columns()]
+
                 self._face_lattice = lattice_from_incidences(
                                     atom_to_facets, facet_to_atoms, ConeFace,
                                     key = id(self))
@@ -2733,6 +2733,51 @@ def facets(self):
         """
         return self.faces(codim=1)
 
+    @cached_method
+    def incidence_matrix(self):
+        r"""
+        Return the incidence matrix.
+
+        .. NOTE::
+
+           The columns correspond to facets/facet normals
+           in the order of :meth:`facet_normals`, the rows
+           correspond to the rays in the order of
+           :meth:`rays`.
+
+        EXAMPLES::
+
+            sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
+            sage: octant.incidence_matrix()
+            [0 1 1]
+            [1 0 1]
+            [1 1 0]
+
+            sage: halfspace = Cone([(1,0,0), (0,1,0), (-1,-1,0), (0,0,1)])
+            sage: halfspace.incidence_matrix()
+            [0]
+            [1]
+            [1]
+            [1]
+            [1]
+
+        TESTS::
+
+            sage: halfspace.incidence_matrix().is_immutable()
+            True
+        """
+        normals = self.facet_normals()
+        incidence_matrix = matrix(ZZ, self.nrays(),
+                                  len(normals), 0)
+
+        for Hindex, normal in enumerate(self.facet_normals()):
+            for Vindex, ray in enumerate(self.rays()):
+                if normal*ray == 0:
+                    incidence_matrix[Vindex, Hindex] = 1
+
+        incidence_matrix.set_immutable()
+        return incidence_matrix
+
     def intersection(self, other):
         r"""
         Compute the intersection of two cones.