In this module, a polyhedron is a convex (possibly unbounded) set in Euclidean space cut out by a finite set of linear inequalities and linear equations. Note that the dimension of the polyhedron can be less than the dimension of the ambient space. There are two complementary representations of the same data:
This describes a polyhedron as the common solution set of a finite number of
The other representation is as the convex hull of vertices (and rays and lines to all for unbounded polyhedra) as generators. The polyhedron is then the Minkowski sum
where
When specifying a polyhedron, you can input a non-minimal set of inequalities/equations or generating vertices/rays/lines. The non-minimal generators are usually called points, non-extremal rays, and non-extremal lines, but for our purposes it is more convenient to always talk about vertices/rays/lines. Sage will remove any superfluous representation objects and always return a minimal representation. For example, \((0,0)\) is a superfluous vertex here:
sage: triangle = Polyhedron(vertices=[(0,2), (-1,0), (1,0), (0,0)])
sage: triangle.vertices()
(A vertex at (-1, 0), A vertex at (1, 0), A vertex at (0, 2))
A polytope is defined as a bounded polyhedron. In this case, the minimal representation is unique and a vertex of the minimal representation is equivalent to a 0-dimensional face of the polytope. This is why one generally does not distinguish vertices and 0-dimensional faces. But for non-bounded polyhedra we have to allow for a more general notion of “vertex” in order to make sense of the Minkowsi sum presentation:
sage: half_plane = Polyhedron(ieqs=[(0,1,0)])
sage: half_plane.Hrepresentation()
(An inequality (1, 0) x + 0 >= 0,)
sage: half_plane.Vrepresentation()
(A line in the direction (0, 1), A ray in the direction (1, 0), A vertex at (0, 0))
Note how we need a point in the above example to anchor the ray and line. But any point on the boundary of the half-plane would serve the purpose just as well. Sage picked the origin here, but this choice is not unique. Similarly, the choice of ray is arbitrary but necessary to generate the half-plane.
Finally, note that while rays and lines generate unbounded edges of the polyhedron they are not in a one-to-one correspondence with them. For example, the infinite strip has two infinite edges (1-faces) but only one generating line:
sage: strip = Polyhedron(vertices=[(1,0),(-1,0)], lines=[(0,1)])
sage: strip.Hrepresentation()
(An inequality (1, 0) x + 1 >= 0, An inequality (-1, 0) x + 1 >= 0)
sage: strip.lines()
(A line in the direction (0, 1),)
sage: strip.faces(1)
(<0,1>, <0,2>)
sage: for face in strip.faces(1):
... print face, '=', face.as_polyhedron().Vrepresentation()
<0,1> = (A line in the direction (0, 1), A vertex at (-1, 0))
<0,2> = (A line in the direction (0, 1), A vertex at (1, 0))
EXAMPLES:
sage: trunc_quadr = Polyhedron(vertices=[[1,0],[0,1]], rays=[[1,0],[0,1]])
sage: trunc_quadr
A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 2 vertices and 2 rays
sage: v = next(trunc_quadr.vertex_generator()) # the first vertex in the internal enumeration
sage: v
A vertex at (0, 1)
sage: v.vector()
(0, 1)
sage: list(v)
[0, 1]
sage: len(v)
2
sage: v[0] + v[1]
1
sage: v.is_vertex()
True
sage: type(v)
<class 'sage.geometry.polyhedron.representation.Vertex'>
sage: type( v() )
<type 'sage.modules.vector_integer_dense.Vector_integer_dense'>
sage: v.polyhedron()
A 2-dimensional polyhedron in ZZ^2 defined as the convex hull of 2 vertices and 2 rays
sage: r = next(trunc_quadr.ray_generator())
sage: r
A ray in the direction (0, 1)
sage: r.vector()
(0, 1)
sage: list( v.neighbors() )
[A ray in the direction (0, 1), A vertex at (1, 0)]
Inequalities \(A \vec{x} + b \geq 0\) (and, similarly, equations) are specified by a list [b, A]:
sage: Polyhedron(ieqs=[(0,1,0),(0,0,1),(1,-1,-1)]).Hrepresentation()
(An inequality (-1, -1) x + 1 >= 0,
An inequality (1, 0) x + 0 >= 0,
An inequality (0, 1) x + 0 >= 0)
See Polyhedron() for a detailed description of all possible ways to construct a polyhedron.
The base ring of the polyhedron can be specified by the base_ring optional keyword argument. If not specified, a suitable common base ring for all coordinates/coefficients will be chosen automatically. Important cases are:
Polyhedra with symmetries often are defined over some algebraic field extension of the rationals. As a simple example, consider the equilateral triangle whose vertex coordinates involve \(\sqrt{3}\). An exact way to work with roots in Sage is the Algebraic Real Field
sage: triangle = Polyhedron([(0,0), (1,0), (1/2, sqrt(3)/2)], base_ring=AA)
sage: triangle.Hrepresentation()
(An inequality (-1, -0.5773502691896258?) x + 1 >= 0,
An inequality (1, -0.5773502691896258?) x + 0 >= 0,
An inequality (0, 1.154700538379252?) x + 0 >= 0)
Without specifying the base_ring, the sqrt(3) would be a symbolic ring element and, therefore, the polyhedron defined over the symbolic ring. This is possible as well, but rather slow:
sage: Polyhedron([(0,0), (1,0), (1/2, sqrt(3)/2)])
A 2-dimensional polyhedron in (Symbolic Ring)^2 defined as the convex
hull of 3 vertices
Even faster than all algebraic real numbers (the field AA) is to take the smallest extension field. For the equilateral triangle, that would be:
sage: K.<sqrt3> = NumberField(x^2-3)
sage: Polyhedron([(0,0), (1,0), (1/2, sqrt3/2)])
A 2-dimensional polyhedron in (Number Field in sqrt3 with defining
polynomial x^2 - 3)^2 defined as the convex hull of 3 vertices
REFERENCES:
Komei Fukuda’s FAQ in Polyhedral Computation
AUTHORS:
- Marshall Hampton: first version, bug fixes, and various improvements, 2008 and 2009
- Arnaud Bergeron: improvements to triangulation and rendering, 2008
- Sebastien Barthelemy: documentation improvements, 2008
- Volker Braun: refactoring, handle non-compact case, 2009 and 2010
- Andrey Novoseltsev: added Hasse_diagram_from_incidences, 2010
- Volker Braun: rewrite to use PPL instead of cddlib, 2011
- Volker Braun: Add support for arbitrary subfields of the reals
Construct a polyhedron object.
You may either define it with vertex/ray/line or inequalities/equations data, but not both. Redundant data will automatically be removed (unless minimize=False), and the complementary representation will be computed.
INPUT:
Some backends support further optional arguments:
OUTPUT:
The polyhedron defined by the input data.
EXAMPLES:
Construct some polyhedra:
sage: square_from_vertices = Polyhedron(vertices = [[1, 1], [1, -1], [-1, 1], [-1, -1]])
sage: square_from_ieqs = Polyhedron(ieqs = [[1, 0, 1], [1, 1, 0], [1, 0, -1], [1, -1, 0]])
sage: list(square_from_ieqs.vertex_generator())
[A vertex at (1, -1),
A vertex at (1, 1),
A vertex at (-1, 1),
A vertex at (-1, -1)]
sage: list(square_from_vertices.inequality_generator())
[An inequality (1, 0) x + 1 >= 0,
An inequality (0, 1) x + 1 >= 0,
An inequality (-1, 0) x + 1 >= 0,
An inequality (0, -1) x + 1 >= 0]
sage: p = Polyhedron(vertices = [[1.1, 2.2], [3.3, 4.4]], base_ring=RDF)
sage: p.n_inequalities()
2
The same polyhedron given in two ways:
sage: p = Polyhedron(ieqs = [[0,1,0,0],[0,0,1,0]])
sage: p.Vrepresentation()
(A line in the direction (0, 0, 1),
A ray in the direction (1, 0, 0),
A ray in the direction (0, 1, 0),
A vertex at (0, 0, 0))
sage: q = Polyhedron(vertices=[[0,0,0]], rays=[[1,0,0],[0,1,0]], lines=[[0,0,1]])
sage: q.Hrepresentation()
(An inequality (1, 0, 0) x + 0 >= 0,
An inequality (0, 1, 0) x + 0 >= 0)
Finally, a more complicated example. Take \(\mathbb{R}_{\geq 0}^6\) with coordinates \(a, b, \dots, f\) and
- The inequality \(e+b \geq c+d\)
- The inequality \(e+c \geq b+d\)
- The equation \(a+b+c+d+e+f = 31\)
sage: positive_coords = Polyhedron(ieqs=[
... [0, 1, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0],
... [0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 0, 1]])
sage: P = Polyhedron(ieqs=positive_coords.inequalities() + (
... [0,0,1,-1,-1,1,0], [0,0,-1,1,-1,1,0]), eqns=[[-31,1,1,1,1,1,1]])
sage: P
A 5-dimensional polyhedron in QQ^6 defined as the convex hull of 7 vertices
sage: P.dim()
5
sage: P.Vrepresentation()
(A vertex at (31, 0, 0, 0, 0, 0), A vertex at (0, 0, 0, 0, 0, 31),
A vertex at (0, 0, 0, 0, 31, 0), A vertex at (0, 0, 31/2, 0, 31/2, 0),
A vertex at (0, 31/2, 31/2, 0, 0, 0), A vertex at (0, 31/2, 0, 0, 31/2, 0),
A vertex at (0, 0, 0, 31/2, 31/2, 0))
Note