fix neighbors of Hrepresentatives

Author: Jonathan Kliem

src/sage/geometry/polyhedron/base.py

@@ -3060,14 +3060,14 @@ def normal_fan(self, direction='inner'):
         r"""
         Return the normal fan of a compact full-dimensional rational polyhedron.
 
-        This returns the inner normal fan of ``self``. For the outer normal fan, 
+        This returns the inner normal fan of ``self``. For the outer normal fan,
         use ``direction='outer'``.
 
         INPUT:
 
-        - ``direction`` -- either ``'inner'`` (default) or ``'outer'``; if 
-          set to ``'inner'``, use the inner normal vectors to span the cones of 
-          the fan, if set to ``'outer'``, use the outer normal vectors. 
+        - ``direction`` -- either ``'inner'`` (default) or ``'outer'``; if
+          set to ``'inner'``, use the inner normal vectors to span the cones of
+          the fan, if set to ``'outer'``, use the outer normal vectors.
 
         OUTPUT:
 
@@ -4388,6 +4388,21 @@ def stack(self, face, position=None):
             ...
             ValueError: the chosen position is too large
 
+            sage: s = polytopes.simplex(7)
+            sage: f = s.faces(3)[0]
+            sage: sf = s.stack(f); sf
+            A 7-dimensional polyhedron in QQ^8 defined as the convex hull of 9 vertices
+            sage: sf.vertices()
+            (A vertex at (-4, -4, -4, -4, 17/4, 17/4, 17/4, 17/4),
+             A vertex at (0, 0, 0, 0, 0, 0, 0, 1),
+             A vertex at (0, 0, 0, 0, 0, 0, 1, 0),
+             A vertex at (0, 0, 0, 0, 0, 1, 0, 0),
+             A vertex at (0, 0, 0, 0, 1, 0, 0, 0),
+             A vertex at (0, 0, 0, 1, 0, 0, 0, 0),
+             A vertex at (0, 0, 1, 0, 0, 0, 0, 0),
+             A vertex at (0, 1, 0, 0, 0, 0, 0, 0),
+             A vertex at (1, 0, 0, 0, 0, 0, 0, 0))
+
         It is possible to stack on unbounded faces::
 
             sage: Q = Polyhedron(vertices=[[0,1],[1,0]],rays=[[1,1]])
@@ -4430,9 +4445,9 @@ def stack(self, face, position=None):
 
         # Taking all facets that contain the face
         if face.dim() == self.dim() - 1:
-            face_star = set([face.ambient_Hrepresentation()[0]])
+            face_star = set([face.ambient_Hrepresentation()[-1]])
         else:
-            face_star = set(facet for facet in self.Hrepresentation()
+            face_star = set(facet for facet in self.Hrepresentation() if facet.is_inequality()
                             if all(facet.contains(x) and not facet.interior_contains(x) for x in face_vertices))
 
         neighboring_facets = set()

src/sage/geometry/polyhedron/representation.py

@@ -449,6 +449,8 @@ def neighbors(self):
         """
         Iterate over the adjacent facets (i.e. inequalities/equations)
 
+        Only defined for inequalities.
+
         EXAMPLES::
 
             sage: p = Polyhedron(ieqs = [[0,0,0,1],[0,0,1,0,],[0,1,0,0],
@@ -459,11 +461,35 @@ def neighbors(self):
             An inequality (0, -1, 0) x + 1 >= 0
             sage: list(a[0])
             [1, 0, -1, 0]
+
+        TESTS:
+
+        Checking that it works for polyhedra with equations::
+
+        sage: P = polytopes.simplex()
+        sage: F1 = P.Hrepresentation()[1]
+        sage: list(F1.neighbors())
+        [An inequality (0, 1, 0, 0) x + 0 >= 0,
+         An inequality (0, 0, 1, 0) x + 0 >= 0,
+         An inequality (0, 0, 0, 1) x + 0 >= 0]
+
+        Does not work for inequalities::
+
+        sage: F0 = P.Hrepresentation()[0]
+        sage: list(F0.neighbors())
+        Traceback (most recent call last):
+        ...
+        AssertionError: must be inequality
         """
+        # The adjacency matrix does not include equations.
+        n_eqs = self.polyhedron().n_equations()
+        assert self.is_inequality(), "must be inequality"
+
         adjacency_matrix = self.polyhedron().facet_adjacency_matrix()
         for x in self.polyhedron().Hrep_generator():
-            if adjacency_matrix[self.index(), x.index()] == 1:
-                yield x
+            if not x.is_equation():
+                if adjacency_matrix[self.index()-n_eqs, x.index()-n_eqs] == 1:
+                    yield x
 
     def adjacent(self):
         """
@@ -777,9 +803,9 @@ def outer_normal(self):
         Return the outer normal vector of ``self``.
 
         OUTPUT:
-        
+
         The normal vector directed away from the interior of the polyhedron.
-        
+
         EXAMPLES::
 
             sage: p = Polyhedron(vertices = [[0,0,0],[1,1,0],[1,2,0]])