Skip to content

Commit

Permalink
Trac #33009: Move is_prism and is_bipyramid to combinatorial polyhedron
Browse files Browse the repository at this point in the history 
This also cleans some code duplication introduced in #27689, as
obtaining the dual of a combinatorial polyhedron is fast:

{{{
sage: %time P = polytopes.hypercube(7)
CPU times: user 31.4 ms, sys: 182 µs, total: 31.6 ms
Wall time: 29.7 ms
sage: %time P.polar()
CPU times: user 25 ms, sys: 7.85 ms, total: 32.8 ms
Wall time: 31.6 ms
A 7-dimensional polyhedron in ZZ^7 defined as the convex hull of 14
vertices
sage: C = P.combinatorial_polyhedron()
sage: %time C.polar()
CPU times: user 57 µs, sys: 0 ns, total: 57 µs
Wall time: 58.4 µs
A 7-dimensional combinatorial polyhedron with 128 facets
}}}

URL: https://trac.sagemath.org/33009
Reported by: gh-kliem
Ticket author(s): Jonathan Kliem
Reviewer(s): Matthias Koeppe
  • Loading branch information
Release Manager committed Feb 12, 2022
2 parents 87011e9 + 80666d4 commit c45f12e
Show file tree
Hide file tree
Showing 2 changed files with 199 additions and 137 deletions.
157 changes: 20 additions & 137 deletions src/sage/geometry/polyhedron/base3.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1276,13 +1276,13 @@ def is_bipyramid(self, certificate=False):
sage: P.is_bipyramid()
True
sage: P.is_bipyramid(certificate=True)
(True, [A vertex at (-1, 0, 0), A vertex at (1, 0, 0)])
(True, [A vertex at (1, 0, 0), A vertex at (-1, 0, 0)])
sage: Q = polytopes.cyclic_polytope(3,7)
sage: Q.is_bipyramid()
False
sage: R = Q.bipyramid()
sage: R.is_bipyramid(certificate=True)
(True, [A vertex at (-1, 3, 13, 63), A vertex at (1, 3, 13, 63)])
(True, [A vertex at (1, 3, 13, 63), A vertex at (-1, 3, 13, 63)])
TESTS::
Expand All @@ -1299,72 +1299,11 @@ def is_bipyramid(self, certificate=False):
Traceback (most recent call last):
...
ValueError: polyhedron has to be compact
ALGORITHM:
Assume all faces of a polyhedron to be given as lists of vertices.
A polytope is a bipyramid with apexes `v`, `w` if and only if for each
proper face `v \in F` there exists a face `G` with
`G \setminus \{w\} = F \setminus \{v\}`
and vice versa (for each proper face
`w \in F` there exists ...).
To check this property it suffices to check for all facets of the polyhedron.
"""
if not self.is_compact():
raise ValueError("polyhedron has to be compact")

from sage.misc.functional import is_odd
n_verts = self.n_vertices()
n_facets = self.n_facets()
if is_odd(n_facets):
if certificate:
return (False, None)
return False

IM = self.incidence_matrix()
if self.n_equations():
# Remove equations from the incidence matrix,
# such that this is the vertex-facet incidences matrix.
I1 = IM.transpose()
I2 = I1[[i for i in range(self.n_Hrepresentation())
if not self.Hrepresentation()[i].is_equation()]]
IM = I2.transpose()

facets_incidences = [set(column.nonzero_positions()) for column in IM.columns()]
verts_incidences = dict()
for i in range(n_verts):
v_i = set(IM.row(i).nonzero_positions())
if len(v_i) == n_facets/2:
verts_incidences[i] = v_i

# Find two vertices ``vert1`` and ``vert2`` such that one of them
# lies on exactly half of the facets, and the other one lies on
# exactly the other half.
from itertools import combinations
for index1, index2 in combinations(verts_incidences, 2):
vert1_incidences = verts_incidences[index1]
vert2_incidences = verts_incidences[index2]
vert1and2 = vert1_incidences.union(vert2_incidences)
if len(vert1and2) == n_facets:
# We have found two candidates for apexes.
# Remove from each facet ``index1`` resp. ``index2``.
test_facets = set(frozenset(facet_inc.difference({index1, index2}))
for facet_inc in facets_incidences)
if len(test_facets) == n_facets/2:
# For each `F` containing `index1` there is
# `G` containing `index2` such that
# `F \setminus \{index1\} = G \setminus \{index2\}
# and vice versa.
if certificate:
V = self.vertices()
return (True, [V[index1], V[index2]])
return True

if certificate:
return (False, None)
return False
return self.combinatorial_polyhedron().is_bipyramid(certificate)

def is_prism(self, certificate=False):
"""
Expand Down Expand Up @@ -1395,36 +1334,36 @@ def is_prism(self, certificate=False):
True
sage: P.is_prism(certificate=True)
(True,
[[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)),
(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),
A vertex at (-1, 1, 1)]])
A vertex at (-1, 1, 1))])
sage: Q = polytopes.cyclic_polytope(3,8)
sage: Q.is_prism()
False
sage: R = Q.prism()
sage: R.is_prism(certificate=True)
(True,
[[A vertex at (1, 6, 36, 216),
A vertex at (1, 0, 0, 0),
A vertex at (1, 7, 49, 343),
A vertex at (1, 5, 25, 125),
A vertex at (1, 1, 1, 1),
A vertex at (1, 2, 4, 8),
A vertex at (1, 4, 16, 64),
A vertex at (1, 3, 9, 27)],
[A vertex at (0, 3, 9, 27),
[(A vertex at (0, 3, 9, 27),
A vertex at (0, 6, 36, 216),
A vertex at (0, 0, 0, 0),
A vertex at (0, 7, 49, 343),
A vertex at (0, 5, 25, 125),
A vertex at (0, 1, 1, 1),
A vertex at (0, 2, 4, 8),
A vertex at (0, 4, 16, 64)]])
A vertex at (0, 4, 16, 64)),
(A vertex at (1, 6, 36, 216),
A vertex at (1, 0, 0, 0),
A vertex at (1, 7, 49, 343),
A vertex at (1, 5, 25, 125),
A vertex at (1, 1, 1, 1),
A vertex at (1, 2, 4, 8),
A vertex at (1, 4, 16, 64),
A vertex at (1, 3, 9, 27))])
TESTS::
Expand All @@ -1441,67 +1380,11 @@ def is_prism(self, certificate=False):
Traceback (most recent call last):
...
NotImplementedError: polyhedron has to be compact
ALGORITHM:
See :meth:`Polyhedron_base.is_bipyramid`.
"""
if not self.is_compact():
raise NotImplementedError("polyhedron has to be compact")

from sage.misc.functional import is_odd
n_verts = self.n_vertices()
n_facets = self.n_facets()
if is_odd(n_verts):
if certificate:
return (False, None)
return False

IM = self.incidence_matrix()
if self.n_equations():
# Remove equations from the incidence matrix,
# such that this is the vertex-facet incidences matrix.
I1 = IM.transpose()
I2 = I1[[i for i in range(self.n_Hrepresentation())
if not self.Hrepresentation()[i].is_equation()]]
IM = I2.transpose()

verts_incidences = [set(row.nonzero_positions()) for row in IM.rows()]
facets_incidences = dict()
for j in range(n_facets):
F_j = set(IM.column(j).nonzero_positions())
if len(F_j) == n_verts/2:
facets_incidences[j] = F_j

# Find two vertices ``facet1`` and ``facet2`` such that one of them
# contains exactly half of the vertices, and the other one contains
# exactly the other half.
from itertools import combinations
for index1, index2 in combinations(facets_incidences, 2):
facet1_incidences = facets_incidences[index1]
facet2_incidences = facets_incidences[index2]
facet1and2 = facet1_incidences.union(facet2_incidences)
if len(facet1and2) == n_verts:
# We have found two candidates for base faces.
# Remove from each vertex ``index1`` resp. ``index2``.
test_verts = set(frozenset(vert_inc.difference({index1, index2}))
for vert_inc in verts_incidences)
if len(test_verts) == n_verts/2:
# For each vertex containing `index1` there is
# another one contained in `index2`
# and vice versa.
# Other than `index1` and `index2` both are contained in
# exactly the same facets.
if certificate:
V = self.vertices()
facet1_vertices = [V[i] for i in facet1_incidences]
facet2_vertices = [V[i] for i in facet2_incidences]
return (True, [facet1_vertices, facet2_vertices])
return True

if certificate:
return (False, None)
return False
return self.combinatorial_polyhedron().is_prism(certificate)

def is_lawrence_polytope(self):
"""
Expand Down
179 changes: 179 additions & 0 deletions src/sage/geometry/polyhedron/combinatorial_polyhedron/base.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2415,6 +2415,185 @@ cdef class CombinatorialPolyhedron(SageObject):
return (False, None)
return False

@cached_method
def is_bipyramid(self, certificate=False):
r"""
Test whether the polytope is a bipyramid over some other polytope.
INPUT:
- ``certificate`` -- boolean (default: ``False``); specifies whether
to return a vertex of the polytope which is the apex of a pyramid,
if found
INPUT:
- ``certificate`` -- boolean (default: ``False``); specifies whether
to return two vertices of the polytope which are the apices of a
bipyramid, if found
OUTPUT:
If ``certificate`` is ``True``, returns a tuple containing:
1. Boolean.
2. ``None`` or a tuple containing:
a. The first apex.
b. The second apex.
If ``certificate`` is ``False`` returns a boolean.
EXAMPLES::
sage: C = polytopes.hypercube(4).combinatorial_polyhedron()
sage: C.is_bipyramid()
False
sage: C.is_bipyramid(certificate=True)
(False, None)
sage: C = polytopes.cross_polytope(4).combinatorial_polyhedron()
sage: C.is_bipyramid()
True
sage: C.is_bipyramid(certificate=True)
(True, [A vertex at (1, 0, 0, 0), A vertex at (-1, 0, 0, 0)])
For unbounded polyhedra, an error is raised::
sage: C = CombinatorialPolyhedron([[0,1], [0,2]], far_face=[1,2], unbounded=True)
sage: C.is_pyramid()
Traceback (most recent call last):
...
ValueError: polyhedron has to be compact
TESTS::
sage: CombinatorialPolyhedron(-1).is_bipyramid()
False
sage: CombinatorialPolyhedron(-1).is_bipyramid(True)
(False, None)
sage: C = polytopes.cross_polytope(1)
sage: C.is_bipyramid()
True
sage: C.is_bipyramid(True)
(True, [A vertex at (1), A vertex at (-1)])
Check that bug analog to :trac:`30292` is avoided::
sage: Polyhedron([[0, 1, 0], [0, 0, 1], [0, -1, -1], [1, 0, 0], [-1, 0, 0]]).is_bipyramid(certificate=True)
(True, [A vertex at (1, 0, 0), A vertex at (-1, 0, 0)])
ALGORITHM:
Assume all faces of a polyhedron to be given as lists of vertices.
A polytope is a bipyramid with apexes `v`, `w` if and only if for each
proper face `v \in F` there exists a face `G` with
`G \setminus \{w\} = F \setminus \{v\}`
and vice versa (for each proper face
`w \in F` there exists ...).
To check this property it suffices to check for all facets of the polyhedron.
"""
if not self.is_compact():
raise ValueError("polyhedron has to be compact")

n_facets = self.n_facets()
if n_facets % 2 or self.dim() < 1:
if certificate:
return (False, None)
return False

facets_incidences = [set(f) for f in self.facets(names=False)]
verts_incidences = dict()
for v in self.face_iter(0):
verts_incidences[v.ambient_V_indices()[0]] = set(v.ambient_H_indices(add_equations=False))

# Find two vertices ``vert1`` and ``vert2`` such that one of them
# lies on exactly half of the facets, and the other one lies on
# exactly the other half.
from itertools import combinations
for index1, index2 in combinations(verts_incidences, 2):
vert1_incidences = verts_incidences[index1]
vert2_incidences = verts_incidences[index2]
vert1and2 = vert1_incidences.union(vert2_incidences)
if len(vert1and2) == n_facets:
# We have found two candidates for apexes.
# Remove from each facet ``index1`` resp. ``index2``.
test_facets = set(frozenset(facet_inc.difference({index1, index2}))
for facet_inc in facets_incidences)
if len(test_facets) == n_facets/2:
# For each `F` containing `index1` there is
# `G` containing `index2` such that
# `F \setminus \{index1\} = G \setminus \{index2\}
# and vice versa.
if certificate:
V = self.vertices()
return (True, [V[index1], V[index2]])
return True

if certificate:
return (False, None)
return False

@cached_method
def is_prism(self, certificate=False):
r"""
Test whether the polytope is a prism of some polytope.
INPUT:
- ``certificate`` -- boolean (default: ``False``); specifies whether
to return two facets of the polytope which are the bases of a prism,
if found
OUTPUT:
If ``certificate`` is ``True``, returns a tuple containing:
1. Boolean.
2. ``None`` or a tuple containing:
a. List of the vertices of the first base facet.
b. List of the vertices of the second base facet.
If ``certificate`` is ``False`` returns a boolean.
TESTS::
sage: CombinatorialPolyhedron(-1).is_prism()
False
sage: CombinatorialPolyhedron(1).is_prism()
False
sage: C = polytopes.cross_polytope(3).prism().combinatorial_polyhedron()
sage: C.is_prism(certificate=True)
(True,
[(A vertex at (0, 0, 1, 0),
A vertex at (0, 1, 0, 0),
A vertex at (0, 0, 0, -1),
A vertex at (0, 0, -1, 0),
A vertex at (0, -1, 0, 0),
A vertex at (0, 0, 0, 1)),
(A vertex at (1, 1, 0, 0),
A vertex at (1, 0, 0, -1),
A vertex at (1, 0, -1, 0),
A vertex at (1, -1, 0, 0),
A vertex at (1, 0, 0, 1),
A vertex at (1, 0, 1, 0))])
sage: C = CombinatorialPolyhedron([[0,1], [0,2]], far_face=[1,2], unbounded=True)
sage: C.is_prism()
Traceback (most recent call last):
...
ValueError: self must be bounded
"""
if not certificate:
return self.dual().is_bipyramid()

val, cert = self.dual().is_bipyramid(True)
if val:
facets = self.facets()
return (True, [facets[cert[0]], facets[cert[1]]])

return (False, None)


def join_of_Vrep(self, *indices):
r"""
Return the smallest face containing all Vrepresentatives indicated by the indices.
Expand Down

0 comments on commit c45f12e

Please sign in to comment.