A Brief Introduction to Polytopes in Sage¶
Author: sarah-marie belcastro <smbelcas@toroidalsnark.net>
If you already know some convex geometry a la Grünbaum or Brøndsted, then you may have itched to get your hands dirty with some polytope calculations. Here is a mini-guide to doing just that.
Basics¶
First, let’s define a polytope as the convex hull of a set of points, i.e. given \(S\) we compute \(P={\rm conv}(S)\):
sage: P1 = Polyhedron(vertices = [[-5,2], [4,4], [3,0], [1,0], [2,-4], [-3,-1], [-5,-3]])
sage: P1
A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 4 vertices
Notice that Sage tells you the dimension of the polytope as well as the dimension of the ambient space.
Of course, you want to know what this object looks like:
sage: P1.plot()
Graphics object consisting of 6 graphics primitives
Even in only 2 dimensions, it’s a pain to figure out what the supporting hyperplanes are. Luckily Sage will take care of that for us.
sage: for q in P1.Hrepresentation():
....: print(q)
An inequality (-4, 1) x + 12 >= 0
An inequality (1, 7) x + 26 >= 0
An inequality (1, 0) x + 5 >= 0
An inequality (2, -9) x + 28 >= 0
That notation is not immediately parseable, because seriously, those do not look like equations of lines (or of halfspaces, which is really what they are).
(-4, 1) x + 12 >= 0 really means \((-4, 1)\cdot\vec{x} + 12 \geq 0\).
So… if you want to define a polytope via inequalities, you have to translate each inequality into a vector. For example, \((-4, 1)\cdot\vec{x} + 12 \geq 0\) becomes (12, -4, 1).
sage: altP1 = Polyhedron(ieqs=[(12, -4, 1), (26, 1, 7),(5,1,0), (28, 2, -9)])
sage: altP1.plot()
Graphics object consisting of 6 graphics primitives
Other information you might want to pull out of Sage about a polytope is the vertex list, which can be done in two ways:
sage: for q in P1.Vrepresentation():
....: print(q)
A vertex at (-5, -3)
A vertex at (-5, 2)
A vertex at (4, 4)
A vertex at (2, -4)
sage: P1.vertices()
(A vertex at (-5, -3), A vertex at (-5, 2), A vertex at (4, 4), A vertex at (2, -4))
Polar duals¶
Surely you want to compute the polar dual:
sage: P1dual = P1.polar()
sage: P1dual
A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 4 vertices
Check it out—we started with an integer-lattice polytope and dualized to a rational-lattice polytope. Let’s look at that.
sage: P1dual.plot()
Graphics object consisting of 6 graphics primitives
sage: P1.plot() + P1dual.plot()
Graphics object consisting of 12 graphics primitives
Oh, yeah, unless the polytope is unit-sphere-sized, the dual will be a very different size. Let’s rescale.
sage: ((1/4)*P1).plot() + (4*P1dual).plot()
Graphics object consisting of 12 graphics primitives
If you think that looks a little bit shady, you’re correct. Here is an example that makes the issue a bit clearer.
sage: P2 = Polyhedron(vertices = [[-5,0], [-1,1], [-2,0], [1,0], [-2,-1], [-3,-1], [-5,-1]])
sage: P2
A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 5 vertices
sage: P2dual = P2.polar(); P2dual
A 2-dimensional polyhedron in QQ^2 defined as the convex hull of 5 vertices
sage: P2.plot() + P2dual.plot()
Graphics object consisting of 14 graphics primitives
That is clearly not computing what we think of as the polar dual. But look at this…
sage: P2.plot() + (-1*P2dual).plot()
Graphics object consisting of 14 graphics primitives
Here is what’s going on.
If a polytope P is in \(\ZZ\), then…
- …the dual is inverted in some way, which is vertically for polygons.
- …the dual is taken of P itself.
- …if the origin is not in P, then an error is returned.
However, if a polytope is not in \(\ZZ\), for example if it’s in \(\QQ\) or
RDF, then…
(1’) …the dual is not inverted.
(2’) …the dual is taken of P-translated-so-barycenter-is-at-origin.
Keep all of this in mind as you take polar duals.
Polytope Constructions¶
Minkowski sums! Now with two syntaxes!
sage: P1+P2
A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 8 vertices
sage: P1.minkowski_sum(P2)
A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 8 vertices
Okay, fine. We should have some 3-dimensional examples, at least. (Note that in order to display polytopes effectively you’ll need visualization software such as Javaview and Jmol installed.)
sage: P3 = Polyhedron(vertices=[(0,0,0), (0,0,1/2), (0,1/2,0), (1/2,0,0), (3/4,1/5,3/2)]); P3
A 3-dimensional polyhedron in QQ^3 defined as the convex hull of 5 vertices
sage: P4 = Polyhedron(vertices=[(-1,1,0),(1,1,0),(-1,0,1), (1,0,1),(0,-1,1),(0,1,1)]); P4
A 3-dimensional polyhedron in ZZ^3 defined as the convex hull of 6 vertices
sage: P3.plot() + P4.plot()
Graphics3d Object
sage: (P3+P4).plot()
Graphics3d Object
We can also find the intersection of two polytopes… and this too has two syntaxes!
sage: int12 = P1.intersection(P2*.5); int12.plot()
Graphics object consisting of 7 graphics primitives
sage: int34 = P3 & P4; int34.plot()
Graphics3d Object
Should one wish to translate, one can.
sage: transP2 = P2.translation([2,1])
sage: P2.plot() + transP2.plot()
Graphics object consisting of 14 graphics primitives
Then of course we can take prisms, pyramids, and bipyramids of polytopes…
sage: P2.prism().plot()
Graphics3d Object
sage: P1.pyramid().plot()
Graphics3d Object
sage: P2dual.bipyramid().plot()
Graphics3d Object
Okay, fine. Yes, Sage has some kinds of polytopes built in.
If you type polytopes. and then press TAB after the period, you’ll get a
list of pre-built polytopes.
sage: P5 = polytopes.hypercube(5)
sage: P6 = polytopes.cross_polytope(3)
sage: P7 = polytopes.simplex(7)
Let’s look at a 4-dimensional polytope.
sage: P8 = polytopes.hypercube(4)
sage: P8.plot()
Graphics3d Object
We can see it from a different perspective:
sage: P8.schlegel_projection([2,5,11,17]).plot()
Graphics3d Object
Queries to polytopes¶
Once you’ve constructed some polytope, you can ask Sage questions about it.
sage: P1.contains([1,0])
True
sage: P1.interior_contains([3,0])
False
sage: P3.contains([1,0,0])
False
Face information can be useful.
sage: int34.f_vector()
(1, 8, 12, 6, 1)
Well, geometric information might be more helpful… Here we are told which of the vertices form each 2-face:
sage: [f.ambient_V_indices() for f in int34.faces(2)]
[(1, 3, 4), (0, 1, 3, 5), (0, 1, 2, 4, 6), (2, 3, 4, 5, 7), (2, 6, 7), (0, 5, 6, 7)]
Yeah, that isn’t so useful as it is. Let’s figure out the vertex and hyperplane representations of the first face in the list.
sage: first2faceofint34 = P3.faces(2)[0]
sage: first2faceofint34.ambient_Hrepresentation(); first2faceofint34.vertices()
(An inequality (1, 0, 0) x + 0 >= 0,)
(A vertex at (0, 0, 0), A vertex at (0, 0, 1/2), A vertex at (0, 1/2, 0))
If you want more… Base class for polyhedra is the first place you want to go.