from typing import Optional, List, Union, Tuple, Sequence, Literal
from .solve_lp import solve_lp
import cvxpy as cp
from matplotlib.axes import Axes
from pyvista import Plotter
import numpy as np
from numpy.typing import NDArray
from scipy.spatial import ConvexHull
from scipy.linalg import null_space
from .convexset import ConvexSet, AffineSet
from mpt4py.tolerances import tolerances
from .visualization import PlotterProtocol
from mpt4py import exceptions as mpterr
from .convexhull import get_convex_hull_solver
from mpt4py.base import Scalar, Vector, Matrix, HData, VData
from .interface_cddlib import cdd_fourier_elimination, cdd_redundancy_elimination
[docs]
class Polyhedron(ConvexSet):
"""
Represents a polyhedron - a set defined by the intersection of a finite number of inequalities.
"""
def __init__(self, H: Optional[HData] = None, V: Optional[VData] = None):
"""
A x <= b, Ae x = be : Inequality form
x = V' * v + R' * w, v >=0, w >= 0, <1, v> = 1 : Vertex form
"""
ConvexSet.__init__(self)
self._Hrep = H
self._Vrep = V
self._chebyshev_center = np.nan * np.zeros((self.dim,))
self._chebyshev_radius = np.nan
self._affine_dim = np.nan
if not (self.hasHrep or self.hasVrep):
raise ValueError("Either H-representation or V-representation must be provided")
def __hash__(self) -> int:
"""
Compute a hash value for the Polyhedron object.
The hash is based on the H-representation if available; otherwise, the V-representation is used.
"""
if self.hasHrep:
return self.Hrep.__hash__()
elif self.hasVrep:
return self.Vrep.__hash__()
else:
raise ValueError("Polyhedron must have either H-representation or V-representation to compute hash.")
[docs]
@classmethod
def from_Hrep(cls, *args, **kwargs) -> "Polyhedron":
"""
Create a polyhedron from its H-representation.
This method initializes a polyhedron object using the H-representation
data provided as arguments, which is encapsulated in an :class:`HData` object.
Args:
*args: Positional arguments to initialize the H-representation.
**kwargs: Keyword arguments to initialize the H-representation.
Returns:
An instance of the :class:`Polyhedron` class initialized with
the given H-representation data.
.. Seealso::
- :func:`from_Vrep` - Create a polyhedron instance from its V-representation.
"""
hrep = HData(*args, **kwargs)
return cls(H=hrep)
[docs]
@classmethod
def from_Vrep(cls, *args, **kwargs) -> "Polyhedron":
"""
Create a polyhedron from its V-representation.
This method initializes a polyhedron object using the V-representation
data provided as arguments, which is encapsulated in an :class:`VData` object.
Args:
*args: Positional arguments to initialize the V-representation.
**kwargs: Keyword arguments to initialize the V-representation.
Returns:
An instance of the :class:`Polyhedron` class initialized with
the given V-representation data.
.. Seealso::
- :func:`from_Hrep` - Create a polyhedron instance from its H-representation.
"""
vrep = VData(*args, **kwargs)
return cls(V=vrep)
[docs]
@classmethod
def empty_set(cls, dim: int) -> "Polyhedron":
r"""Create an empty polyhedron of given dimension which has a H-rep :math:`P = \{x \mid \mathbf{0}^{\top} x \leq -1\}`
and a V-rep (no vertices or rays).
"""
Hrep = HData(A=np.zeros((1, dim)), b=np.array([-1]), is_minimal=True)
Vrep = VData(V=np.empty((0, dim)), R=np.empty((0, dim)), is_minimal=True)
return cls(H=Hrep, V=Vrep)
[docs]
@classmethod
def singleton(cls, point: Vector) -> "Polyhedron":
r"""Create a polyhedron that contains only a single point. The H-rep is :math:`P = \{x \mid x = v\}`
where :math:`v` is the given point. The V-rep contains the point as the only vertex and no rays."""
point = np.asarray(point).reshape((-1,))
dim = point.shape[0]
# compared to creating a singleton using from_Vrep, this method also provides a Hrep,
# which can accelerate some polyhedron operations that will be potentially performed
Hrep = HData(Ae=np.eye(dim), be=point, is_minimal=True)
Vrep = VData(V=point.reshape((1, -1)), is_minimal=True)
return cls(H=Hrep, V=Vrep)
[docs]
@classmethod
def from_bounds(cls, lb: Vector, ub: Vector) -> "Polyhedron":
r"""Create a polyhedron from upper and lower bounds.
If the upper bounds contain ``numpy.inf`` or the lower bounds
contain ``-numpy.inf``, the corresponding dimensions are treated as **unbounded**.
The resulting polyhedron has both H-representation and V-representation. The H-rep is
:math:`P = \{x \mid x_{lb} \leq x \leq x_{ub}\}` and the minimal V-rep is computed accordingly.
"""
lb = np.asarray(lb, dtype=float).reshape(-1)
ub = np.asarray(ub, dtype=float).reshape(-1)
if lb.size != ub.size:
raise ValueError("lower bounds and upper bounds must have the same dimension.")
n = lb.size # dimension
if np.isnan(lb).any():
raise ValueError("Lower bounds cannot contain NaN.")
if np.isnan(ub).any():
raise ValueError("Upper bounds cannot contain NaN.")
if np.isneginf(ub).any():
raise ValueError("Upper bounds cannot contain -inf.")
if np.isposinf(lb).any():
raise ValueError("Lower bounds cannot contain +inf.")
# Infeasibility check
infeasible = (lb - ub) >= tolerances['zero']
if np.any(infeasible):
return Polyhedron.empty_set(n)
## Compute H-rep
is_finite_lb = np.isfinite(lb)
is_finite_ub = np.isfinite(ub)
is_eq = is_finite_lb & is_finite_ub & (np.abs(lb - ub) <= tolerances['zero'])
idx_eq = np.where(is_eq)[0]
idx_both = np.where(is_finite_lb & is_finite_ub & ~is_eq)[0] # two-sided finite
idx_lb_only = np.where(is_finite_lb & ~is_finite_ub)[0] # lower-bounded, +inf above
idx_ub_only = np.where(~is_finite_lb & is_finite_ub)[0] # upper-bounded, -inf below
idx_free = np.where(~is_finite_lb & ~is_finite_ub)[0] # completely free
A_rows = []
b_vals = []
Ae_rows = []
be_vals = []
# equalities x_i = v
for i in idx_eq:
ei = np.zeros(n); ei[i] = 1.0
Ae_rows.append(ei)
be_vals.append(0.5 * (lb[i] + ub[i])) # average for numerical safety
# finite upper bounds: e_i^T x <= ub_i
for i in np.where(is_finite_ub & ~is_eq)[0]:
ei = np.zeros(n); ei[i] = 1.0
A_rows.append(ei)
b_vals.append(ub[i])
# finite lower bounds: -e_i^T x <= -lb_i (i.e., x_i >= lb_i)
for i in np.where(is_finite_lb & ~is_eq)[0]:
ei = np.zeros(n); ei[i] = -1.0
A_rows.append(ei)
b_vals.append(-lb[i])
A = np.vstack(A_rows) if A_rows else np.zeros((0, n))
b = np.asarray(b_vals, float) if b_vals else np.zeros((0,))
Ae = np.vstack(Ae_rows) if Ae_rows else np.zeros((0, n))
be = np.asarray(be_vals, float) if be_vals else np.zeros((0,))
Hrep = HData(A=A, b=b, Ae=Ae, be=be, is_minimal=True)
## Compute V-rep
# Rays: where bounds are infinite
rays = []
# +e_i for lower-bounded only (ub = +inf)
for i in idx_lb_only:
ri = np.zeros(n); ri[i] = 1.0
rays.append(ri)
# -e_i for upper-bounded only (lb = -inf)
for i in idx_ub_only:
ri = np.zeros(n); ri[i] = -1.0
rays.append(ri)
# ±e_i (lineality) for completely free
for i in idx_free:
rpos = np.zeros(n); rpos[i] = 1.0
rneg = np.zeros(n); rneg[i] = -1.0
rays.append(rpos)
rays.append(rneg)
R = np.vstack(rays) if rays else np.zeros((0, n))
# Vertices:
# - If any completely free dimension exists (±inf), the set has lineality -> no vertices.
# - Otherwise, take Cartesian product over two-sided finite dims {lb_i, ub_i},
# and fix equality dims at their value, and one-sided dims at their finite bound.
if idx_free.size > 0:
V = np.zeros((0, n))
else:
# Start from a single prototype
base = np.zeros(n)
# equality dims fixed
for i in idx_eq:
base[i] = 0.5 * (lb[i] + ub[i])
# one-sided dims fixed at the finite bound (this is where an extreme point can occur)
for i in idx_lb_only:
base[i] = lb[i] # only lower finite -> fix to lb for a vertex
for i in idx_ub_only:
base[i] = ub[i] # only upper finite -> fix to ub
# For two-sided finite dims, enumerate corners
if idx_both.size == 0:
V = base.reshape(1, -1) # one vertex (or still unbounded via rays)
else:
# 2^k vertices where k = number of two-sided finite dims
choices = []
for i in idx_both:
choices.append(np.array([lb[i], ub[i]], dtype=float))
# Cartesian product
grids = np.meshgrid(*choices, indexing='ij')
num_vertices = 1
for c in choices:
num_vertices *= c.size
V = np.tile(base, (num_vertices, 1))
# fill the varying coords
flat_coords = [g.reshape(-1) for g in grids]
for j, i in enumerate(idx_both):
V[:, i] = flat_coords[j]
Vrep = VData(V=V, R=R, is_minimal=True)
return Polyhedron(H=Hrep, V=Vrep)
@property
def dim(self) -> int:
r"""
Return the polyhedron's dimension of representation.
.. Seealso::
- :func:`affine_dimension` - Compute the affine dimension of the polyhedron.
- :func:`codimension` - Compute the codimension of the polyhedron.
"""
if self._Hrep is not None:
return self._Hrep.dim
if self._Vrep is not None:
return self._Vrep.dim
raise ValueError("Either H-representation or V-representation must be provided")
@property
def hasHrep(self) -> bool:
"""
Tell if the polyhedron has a H-representation or not.
"""
return self._Hrep is not None
@property
def hasVrep(self) -> bool:
"""
Tell if the polyhedron has a V-representation or not.
"""
return self._Vrep is not None
def _getHrep(self, method: Optional[str] = None) -> HData:
"""
Compute the H-representation of the polyhedron from the V-representation
"""
if self._Hrep is None:
solver = get_convex_hull_solver(method)
if self._Vrep is None:
raise ValueError("V-representation must be provided in order to compute the H-representation")
self._Hrep = solver.convexhull(self._Vrep)
return self._Hrep
def _getVrep(self, method: Optional[str] = None) -> VData:
"""
Compute the V-representation of the polyhedron from the H-representation
"""
if self._Vrep is None:
solver = get_convex_hull_solver(method)
if self._Hrep is None:
raise ValueError("H-representation must be provided in order to compute the V-representation")
self._Vrep = solver.enumerate_vertices(self._Hrep)
return self._Vrep
@property
def A(self) -> Matrix:
r""" Return the :math:`A` matrix in H-representation. Will be computed if needed."""
return self._getHrep().A
@property
def b(self) -> Vector:
r""" Return the :math:`b` vector in H-representation. Will be computed if needed."""
return self._getHrep().b
@property
def H(self) -> Matrix:
r""" Return the contatenated :math:`[A \mid b]` matrix in H-representation. Will be computed if needed."""
# Return [A b]
return np.hstack([self.A, self.b.reshape((-1,1))])
@property
def He(self) -> Matrix:
r""" Return the contatenated :math:`[A_e \mid b_e]` matrix in H-representation. Will be computed if needed."""
# Return [Ae be]
return np.hstack([self.Ae, self.be.reshape((-1,1))])
@property
def Ae(self) -> Matrix:
r""" Return the :math:`A_e` matrix in H-representation. Will be computed if needed."""
return self._getHrep().Ae
@property
def be(self) -> Vector:
r""" Return the :math:`b_e` vector in H-representation. Will be computed if needed."""
return self._getHrep().be
@property
def V(self) -> Matrix:
r""" Return the :math:`V` matrix in V-representation. Will be computed if needed."""
return self._getVrep().V
@property
def R(self) -> Matrix:
r""" Return the :math:`R` matrix in V-representation. Will be computed if needed."""
return self._getVrep().R
[docs]
def minHrep(self, inplace: bool = True, use_cdd: bool = False) -> HData:
r"""
Compute the minimal H-representation of the polyhedron.
Args:
inplace: If True, the current object is modified.
use_cdd: If True, the `pycddlib` is used for the redundancy elimination. Otherwise the built-in implementation is used.
Returns:
The minimal H-representation of the polyhedron.
Note:
- If an H-representation is already known, then this function performs redundancy elimination. Otherwise, facet enumeration will be performed first.
- If the polyhedron is empty, no redundancy elimination will be performed.
.. admonition:: Math
The redundancy elimination of H-representation performs the following steps.
First, attempt to identify some irredundant inequalities using ray shooting heuristics if the polyhedron is **full-dimensional**.
Draw a number of rays from a strict interior point along random directions. For each ray, the inequality it intersects first must be irredundant.
Second, identify redundant inequalities using Linear Programming. For each inequality :math:`a_j x \leq b_j`, solve the following LP:
.. math::
\begin{align*}
\max_x &\quad a_j x \\
\text{s.t.} &\quad a_i x \leq b_i, \quad \forall i\neq j \\
&\quad a_j x \leq b_j + 1, \\
&\quad A_e x = b_e. \\
\end{align*}
This LP is guaranteed to be bounded. The inequality :math:`a_j x \leq b_j` is **irredundant** if the LP's optimal objective value is strictly greater than :math:`b_j` and **redundant** otherwise.
"""
_ = self._getHrep()
if self._Hrep.is_minimal:
return self._Hrep
self._Hrep.preprocess()
if use_cdd:
minhrep, _, _ = cdd_redundancy_elimination(self._Hrep)
self._Hrep = minhrep if inplace else self._Hrep
return minhrep
if self.is_empty:
return self._Hrep
# Removing constraints with a zero normal is already done in the HData preprocessing
# # Remove constraints with a zero normal
# mask = np.linalg.norm(A, axis=1) > tolerances['zero_normal']
# A, b = A[mask], b[mask]
ind_ineq_irredundant = []
ind_ineq_redundant = []
ind_ineq_implicit_linearity = [] # TODO: the indices of inequalities that define implicit equalities
ind_ineq_totest = [*range(self.A.shape[0])]
# ----------------------------------------------------------------------
# HEURISTICS TO IDENTIFY SOME REDUNDANT/IRREDUNDANT INEQUALITIES FIRST
# Detect some irredundant inequalities using ray shooting heuristics, when the polyhedron is full-dimensional
use_rayshooting_heuristics = True
if use_rayshooting_heuristics and self.is_full_dim:
x0, strict_interior = self.interior_point()
if strict_interior:
if np.any(np.isinf(x0)) or np.any(np.isnan(x0)):
raise RuntimeError("The interior point contains inf or nan... this should not happen")
if np.any(self.A @ x0 >= self.b - tolerances['zero']) and strict_interior:
raise RuntimeError("The interior point is not strictly feasible...")
num_ray_shooting = 2 * self.A.shape[0] # TODO: number of shooting rays can be tuned
directions = np.random.randn(num_ray_shooting, self.dim) # generate all random directions
directions /= (np.linalg.norm(directions, axis=1, keepdims=True) + 1E-10)
denominators = self.A @ directions.T
# create mask for valid directions (non-zero denominators)
valid_mask = (np.linalg.norm(denominators, ord=np.inf, axis=0) >= tolerances['zero_normal'])
valid_denominators = denominators[:, valid_mask]
if valid_denominators.size > 0:
numerator = self.b - self.A @ x0 # travel distance calculation. Shape: (num_constraints,)
travel_distance = numerator[:, None] / valid_denominators
travel_distance[valid_denominators <= 0] = np.inf # if no intersection along positive direction, set distance to inf
min_indices = np.argmin(travel_distance, axis=0)
min_values = travel_distance[min_indices, np.arange(travel_distance.shape[1])]
valid_min_indices = min_indices[min_values < np.inf]
ind_ineq_irredundant.extend(valid_min_indices.tolist())
ind_ineq_totest = list(set(ind_ineq_totest).difference(set(ind_ineq_irredundant)))
# Remove duplicates (kept your original unique logic)
ind_ineq_irredundant = np.unique(ind_ineq_irredundant).tolist()
# ----------------------------------------------------------------------
# DETECT REDUNDANT/IRREDUNDANT INEQUALITIES BY SOLVING LINEAR PROGRAMS
# Remove redundant constraints
# We solve the problem
# max c'x s.t. Ax <= b and Ae x = be
# where c is the normal of the constraint under test
USE_CVXPY = False
if USE_CVXPY:
x = cp.Variable(self.dim)
c = cp.Parameter(self.dim)
slack = cp.Parameter(self.A.shape[0])
con_eq: List[cp.Constraint] = [self.Ae @ x == self.be] if self.Ae.size > 0 else []
con_ineq: List[cp.Constraint] = [self.A @ x <= self.b + slack] if self.A.size > 0 else []
# con_box: List[cp.Constraint] = [-1e6 <= x, x <= 1e6] # Add box constraints to avoid numerical issues with unbounded polyhedra
prob = cp.Problem(cp.Maximize(c @ x), con_eq + con_ineq) # + con_box)
while len(ind_ineq_totest) > 0:
test = ind_ineq_totest.pop(0)
c.value = self.A[test]
slack.value = np.zeros(self.A.shape[0])
slack.value[test] = 1. # Relax the constraint under test
# slack.value[ind_ineq_redundant] = np.inf # ! Drop redundant constraints. Probably using np.inf is better than using a large number like 1e8
slack.value[ind_ineq_redundant] = 1e8 # ! Drop redundant constraints. I found using np.inf will cause Clarabel fail, so changed to 1e8. But for Gurobi it's fine.
prob.solve(cp.CLARABEL, verbose=False)
if prob.status == cp.OPTIMAL:
assert isinstance(prob.value, np.floating), "Invalid value returned by the solver"
value = prob.value
if value < self.b[test] + tolerances['redundancy_elimination']:
ind_ineq_redundant.append(test)
else:
ind_ineq_irredundant.append(test)
# use_ray_shooting = False
# if use_ray_shooting and self.is_full_dim:
# # Ray shooting, compute the intersecting point
# x_opt = x.value
# z, _ = self.cheby_center() # an interior point
# # solve equation: y = z+(x_opt-z)*t, a_j*y = b_j
# t_list = []
# for j in ind_totest:
# a_j, b_j = A[j], b[j]
# t = (b_j - a_j @ z) / (a_j @ (x_opt - z))
# if t < 0.:
# t = np.inf
# t_list.append(t)
# if len(t_list) > 0:
# argmin = np.argmin(t_list)
# if t_list[argmin] < 1.:
# ind_irredundant.append(ind_totest[argmin])
# ind_totest.remove(ind_totest[argmin])
elif prob.status == cp.INFEASIBLE:
self._affine_dim = -1 # empty set
return self.Hrep
elif prob.status == cp.UNBOUNDED:
raise ValueError("The polyhedron is unbounded after clipping... this should not happen")
else:
raise RuntimeError(f"The LP solver returned status {prob.status}")
else: # Directly call LP solver
while len(ind_ineq_totest) > 0:
test = ind_ineq_totest.pop(0)
slack = np.zeros(self.A.shape[0])
slack[test] = 1. # Relax the constraint under test
mask = [test] + ind_ineq_totest + ind_ineq_irredundant # only keep untested and irredundant constraints
lp_c = -self.A[test] # solvers we consider below do minimization
lp_A, lp_b = self.Ae, self.be
lp_G, lp_h = self.A[mask, :], (self.b + slack)[mask]
# remove the redundant constraints since giving a large rhs will either cause solver failure or slow down the convergece
# mask = list(set(range(self.A.shape[0])).difference(set(ind_ineq_redundant)))
# mask = list(set().difference(set(ind_ineq_redundant)))
# lp_c = -A[test] # solvers we consider below do minimization
# lp_A, lp_b = self.Ae, self.be
# lp_G, lp_h = self.A[mask, :], (self.b + slack)[mask]
status, x_value, obj_value = solve_lp(c=lp_c, A=lp_A, b=lp_b, G=lp_G, h=lp_h)
if status == 'optimal':
value = -obj_value # remember we negated c
if value < self.b[test] + tolerances['redundancy_elimination']:
ind_ineq_redundant.append(test)
else:
ind_ineq_irredundant.append(test)
elif status == 'infeasible':
self._affine_dim = -1 # empty set
return self.Hrep
elif status == 'unbounded':
raise ValueError("The polyhedron is unbounded after clipping... this should not happen")
else:
raise RuntimeError(f"The LP solver returned status {status}")
# ----------------------------------------------------------------------
# Reduce the equality constraints to minimal form
Ae, be = self.Ae, self.be
if self.Ae.size > 0:
U, s, Vt = np.linalg.svd(self.Ae, full_matrices=False)
rank = np.sum(s > tolerances['zero'])
if rank < self.Ae.shape[0]:
Ae = (s[:, np.newaxis] * Vt).astype(np.float64) # Weird pylance error
be = U.T @ self.be
ind_ineq_redundant.sort()
ind_ineq_irredundant.sort()
H = HData(Ae=Ae, be=be, A=self.A[ind_ineq_irredundant], b=self.b[ind_ineq_irredundant], is_minimal=True)
if inplace:
self._Hrep = H
return H
[docs]
def minVrep(self, inplace: bool = True, use_cdd: bool = False) -> VData:
r"""
Returns the minimal V-representation of the polyhedron.
If inplace is True, then the current object is modified.
.. admonition:: Math
The redundancy elimination of V-representation performs the following steps:
1. Eliminate redundant rays. An extreme ray :math:`r_j` is redundant if it can be written as the non-negative
linear combination of other rays :math:`r_i`, i.e., if the following LP is feasible:
.. math::
\begin{align*}
\min_\mu &\quad 0 \\
\text{s.t.} &\quad r_j = \sum\limits_{i\neq j} \mu_i r_i \\
&\quad \mu_i \geq 0
\end{align*}
2. Eliminate redundant vertices. The method we use is an extension of Section 2.19 in `Fukuda's PolyFAQ
<https://people.inf.ethz.ch/fukudak/Doc_pub/polyfaq220115c.pdf>`_.
A vertex :math:`v_j` is redundant if the following LP is infeasible:
.. math::
\begin{align*}
\min\limits_{z, z_0} &\quad 0 \\
\text{s.t.} &\quad z^\top v_j > z_0 \\
&\quad z^\top v_i \leq z_0 \\
&\quad z^\top r_k < 0
\end{align*}
"""
_ = self._getVrep()
if self.Vrep.is_minimal:
return self._Vrep
self.Vrep.preprocess() # TODO: [bug] with this preprocess function, it cannot do non-inplace operation
if use_cdd:
minvrep, idx_redundant, idx_implicit_linearity = cdd_redundancy_elimination(self.Vrep)
self._Vrep = minvrep if inplace else self._Vrep
return minvrep
ind_rays_irredundant = []
ind_rays_redundant = []
ind_rays_totest = [*range(self.R.shape[0])]
ind_vertices_irredundant = []
ind_vertices_redundant = []
ind_vertices_totest = [*range(self.V.shape[0])]
# ----------------------------------------------------------------------
# HEURISTICS TO IDENTIFY SOME REDUNDANT/IRREDUNDANT VERTICES/RAYS FIRST
if self.hasHrep:
AV_minus_b = self.V @ self.A.T - self.b[np.newaxis, :]
feasible_mask = np.all(AV_minus_b <= tolerances['zero'], axis=1)
redundant_mask = np.all(AV_minus_b <= -tolerances['zero'], axis=1)
if not np.all(feasible_mask):
raise RuntimeError("Inconsistent H-rep and V-rep: some vertices are not feasible.")
ind_vertices_redundant += np.where(redundant_mask)[0].tolist()
ind_vertices_totest = [i for i in ind_vertices_totest if i not in np.where(redundant_mask)[0]]
if self.V.size:
# For a bounded polyhedron, if a vertex is the only maximum / minimum along a direction, it must be irredundant
# (if multiple vertices have the same coordinate at one dimension, then we cannot conclude anything. Just think about a rectange in 2d)
if self.R.size == 0:
for d in range(self.dim):
proj_d = self.V[:, d] # the 1d projection of all vertices along dimension d
max_idx = np.where(np.abs(proj_d - np.max(proj_d)) <= tolerances["zero"])[0]
min_idx = np.where(np.abs(proj_d - np.min(proj_d)) <= tolerances["zero"])[0]
if len(max_idx) == 1 and max_idx[0] in ind_vertices_totest: # the latter condition is to avoid adding duplicates or removing indices already removed
ind_vertices_irredundant.append(max_idx[0])
ind_vertices_totest.remove(max_idx[0])
if len(min_idx) == 1 and min_idx[0] in ind_vertices_totest: # the latter condition is to avoid adding duplicates or removing indices already removed
ind_vertices_irredundant.append(min_idx[0])
ind_vertices_totest.remove(min_idx[0])
# Project to lower dimensions to compute convex hull and find some redundant vertices
# vertices that are redundant in the projection are also redundant in the original space
if self.dim >= 3:
from itertools import combinations
for lower_dims in combinations(range(self.dim), 2):
projected_V = self.V[np.ix_(ind_vertices_totest, lower_dims)]
def _find_vertex_indices(V1: np.ndarray, V2: np.ndarray):
# we assume V2 is a subset of V1, and we want to find the indices of V2 in V1
indices = []
for v in V1:
dists = np.linalg.norm(V2 - v, axis=1)
idx = np.where(dists <= tolerances['zero'])[0]
if len(idx) == 0:
raise ValueError(f"Vertex {v} not found in V2 within tolerance {tolerances['zero']}")
# If multiple, take the first one
indices.append(idx[0])
return np.array(indices, dtype=int)
try:
projected_hull = ConvexHull(projected_V)
ind_irredundant_vertices_in_projected_hull = _find_vertex_indices(projected_hull.vertices, projected_V)
ind_redundant_vertices_in_projected_hull = list(set(range(projected_V.shape[0])).difference(set(ind_irredundant_vertices_in_projected_hull.tolist())))
# If a polyhedron has both rays and vertices, the vertices that are redundant in the convex hull
# of vertices are certainly redundant in the original polyhedron. But we cannot say the vertices
# that are irredundant in the convex hull of vertices are irredundant in the original polyhedron.
ind_vertices_redundant.extend(ind_redundant_vertices_in_projected_hull)
ind_vertices_totest.remove(ind for ind in ind_vertices_redundant)
if self.R.size == 0:
ind_vertices_irredundant.extend(ind_irredundant_vertices_in_projected_hull)
ind_vertices_totest.remove(ind for ind in ind_vertices_irredundant)
except Exception as e:
# If the convex hull computation fails, skip this projection
# warnings.warn(f"Convex hull computation failed for projection onto dimensions {lower_dims}: {e}")
continue
# END OF HEURISTICS
# ----------------------------------------------------------------------
# ----------------------------------------------------------------------
# REDUNDANCY ELIMINATION BY SOLVING LINEAR PROGRAMS
# Repeated rows in V and R are already removed in __init__ method of VData, no need to handle it again
USE_CVXPY = False # TODO: make it an argument?
# Remove redundant rays
# Solve the LP to check if ray i is redundant
# min 0
# s.t. R_other' * mu = r_i
# mu >= 0
# where R_other is R without the i-th row
# If feasible, ray i is redundant; otherwise, it is irredundant
if len(ind_rays_totest) > 0:
num_rays = len(ind_rays_totest)
if USE_CVXPY:
# Define reusable LP structure using parameters
target_ray = cp.Parameter(self.dim) # ray to test (r_i)
other_rays = cp.Parameter((num_rays-1, self.dim)) # other rays (excluding r_i)
mu = cp.Variable(num_rays - 1) # coefficients for the combination of other rays
# LP: Check if target_ray = other_rays * mu, with mu >= 0
constraints = [other_rays.T @ mu == target_ray, mu >= 0]
prob = cp.Problem(cp.Minimize(cp.Constant(0)), constraints) # Feasibility problem
for i in range(num_rays):
# Extract the ray to test (r_i) and the remaining rays
r_i = self.R[i, :]
R_other = np.delete(self.R, i, axis=0) # remove r_i # TODO: check this
# Update parameters and solve the LP
target_ray.value = r_i
other_rays.value = R_other
# prob.solve(solver=cp.CLARABEL, verbose=False)
prob.solve(solver=cp.PIQP, verbose=False)
# If feasible, ray is redundant; otherwise, keep it
if prob.status in [cp.INFEASIBLE, cp.UNBOUNDED]:
ind_rays_irredundant.append(i)
elif prob.status in [cp.OPTIMAL]:
ind_rays_redundant.append(i)
else:
raise RuntimeError(f"Numerical error. LP solver returns {prob.status}")
ind_rays_totest.remove(i)
else: # not USE_CVXPY, directly call LP solvers
while len(ind_rays_totest) > 0:
test = ind_rays_totest.pop(0)
r_test = self.R[test, :]
R_other = self.R[ind_rays_totest + ind_rays_irredundant, :] # remove r_test and redundant rays
m = R_other.shape[0]
status, x_opt, obj_opt = solve_lp(
c=np.zeros(m),
A=R_other.T,
b=r_test,
G=-np.eye(m),
h=np.zeros((m,)),
)
if status == 'optimal':
ind_rays_redundant.append(test)
elif status in ['infeasible', 'unbounded']:
ind_rays_irredundant.append(test)
else:
raise RuntimeError(f"Numerical error. LP solver returns {status}")
assert len(ind_rays_totest) == 0, "All rays should have been tested"
irredundant_rays = self.R[ind_rays_irredundant, :]
# Remove redundant vertices
# The implementation is based on Sec. 2.19 in Fukuda's Polyhedral Computation FAQ
if len(ind_vertices_totest) > 0:
num_vertices = len(ind_vertices_totest)
if USE_CVXPY:
# Define reusable LP structure using parameters
target_vertex = cp.Parameter(self.dim) # ray to test (r_i)
other_vertices = cp.Parameter((num_vertices - 1, self.dim)) # other rays (excluding r_i)
z0 = cp.Variable()
z = cp.Variable(self.dim)
objective = cp.vdot(target_vertex, z) - z0
constraints = [other_vertices @ z - z0 <= 0 ] # + tolerances['redundancy_elimination']]
constraints += [cp.vdot(target_vertex, z) - z0 <= 1]
if len(ind_rays_irredundant) > 0:
constraints += [irredundant_rays @ z <= -tolerances['redundancy_elimination']]
prob = cp.Problem(cp.Maximize(objective), constraints)
tested = []
for test in ind_vertices_totest:
# Extract the vertex to test and the remaining vertices
v_test = self.V[test]
V_other = self.V[list(set(ind_vertices_totest) - {test}), :]
# Update parameters and solve the LP
target_vertex.value = v_test
other_vertices.value = V_other
# prob.solve(solver=cp.CLARABEL, verbose=False)
prob.solve(solver=cp.PIQP, verbose=False)
if not prob.status == cp.OPTIMAL:
raise RuntimeError("The LP is not solved to optimality.")
if objective.value > tolerances['redundancy_elimination']:
ind_vertices_irredundant.append(test)
else:
ind_vertices_redundant.append(test)
tested.append(test)
ind_vertices_totest = list(set(ind_vertices_totest).difference(set(tested)))
else: # not USE_CVXPY, directly call LP solvers
while len(ind_vertices_totest) > 0:
test = ind_vertices_totest.pop(0)
v_test = self.V[test, :]
V_other = self.V[ind_vertices_totest + ind_vertices_irredundant, :]
# min 0
# s.t. vi = V_other' * lambda + R' * mu
# lambda >= 0, mu >= 0
# 1' * lambda = 1
lp_A = np.vstack((
np.hstack((V_other.T, self.R[ind_rays_irredundant].T)) if len(ind_rays_irredundant) > 0 else V_other.T,
np.hstack((np.ones((1, V_other.shape[0])), np.zeros((1, len(ind_rays_irredundant))))),
))
lp_b = np.concatenate((v_test, np.array([1.])))
lp_G = -np.eye(V_other.shape[0] + len(ind_rays_irredundant))
lp_h = np.zeros((V_other.shape[0] + len(ind_rays_irredundant), ))
status, _, _ = solve_lp(
c = np.zeros(V_other.shape[0] + len(ind_rays_irredundant)),
A = lp_A, b = lp_b, G = lp_G, h = lp_h,
)
if status == 'optimal':
ind_vertices_redundant.append(test)
elif status == 'infeasible':
ind_vertices_irredundant.append(test)
else:
raise RuntimeError(f"LP solver returns {status}")
assert len(ind_vertices_totest) == 0, "All vertices should have been tested"
irredundant_vertices = self.V[ind_vertices_irredundant, :]
min_vrep = VData(V=irredundant_vertices, R=irredundant_rays, is_minimal=True)
if inplace:
self._Vrep = min_vrep
return min_vrep
@property
def Vrep(self) -> VData:
"""
Returns the V-representation of the polyhedron. Computes it if needed.
"""
return self._getVrep()
@property
def Hrep(self) -> HData:
"""
Returns the H-representation of the polyhedron. Computes it if needed.
"""
return self._getHrep()
@property
def incidence_map(self) -> dict:
r"""
Compute the incidence map of the polyhedron.
The incidence map describes containment of vertices in facets.
.. Note::
This function computes both irredundant H-rep and V-rep and can be time consuming.
.. warning::
This function is **not implemented yet**. Calling it will raise
``NotImplementedError``.
"""
raise NotImplementedError
Hrep = self.minHrep(inplace=False)
Vrep = self.minVrep(inplace=False)
import scipy.sparse as sp
map = {
'V': Vrep.V,
'R': Vrep.R,
'H': Hrep.H,
'He': Hrep.He,
'incVH': [],
'incRH': []}
# val = self.A @ self.V.T - self.b.reshape((-1,1))
# return abs(val) < tolerances['zero']
# Returns a matrix I in [bool]^{num_inequ x num_points}.
# I[i,j] = 1 if point j is incident on facet i, 0 otherwise.
def _cvxrep(self) -> Tuple[cp.Variable, List[cp.Constraint]]:
"""
Returns a CVXPY representation of the set.
Returns a variable and a set of constraints such that
x in Set <=> x satisfies constraints
"""
x = cp.Variable(self.dim)
if self.hasHrep:
constraints = [self.A @ x <= self.b]
if self.Ae.size > 0:
constraints += [self.Ae @ x == self.be]
return x, constraints
if self.hasVrep:
constraints = []
var = np.zeros((self.dim, ))
if self.V.shape[0] > 0:
v = cp.Variable(self.V.shape[0], nonneg=True)
var = var + self.V.T @ v
constraints += [cp.sum(v) == 1]
if self.R.shape[0] > 0:
w = cp.Variable(self.R.shape[0], nonneg=True)
var = var + self.R.T @ w
constraints += [x == var]
return x, constraints
raise ValueError("Invalid representation - didn't find either an H-representation or a V-representation")
def __contains__(self, inner: Union[Vector, ConvexSet]) -> bool:
"""
Return True if the given point or ConvexSet is contained in this polyhedron and False otherwise
"""
return self.contains(inner)
[docs]
def contains(self, inner: Union[Vector, ConvexSet]) -> bool:
r"""Determine if the given point or convex set is a subset of this polyhedron or not.
Note:
- The H-representation of the self polyhedron is required to test if one set is contained in another. Will be computed if needed.
- If the inner set is unbounded and not a polyhedron, then this will likely raise an error.
.. admonition:: Math
For two homogenized polyhedra :math:`P = \{x \mid Ax\leq 0, A_e x=0 \}` and :math:`Q=\{x\mid Bx\leq 0, B_e x=0 \}`, to test if :math:`Q\subseteq P`:
- If :math:`Q` has V-rep, just test that all vertices and rays of :math:`Q` are contained in :math:`P`.
- If :math:`Q` only has H-rep, use Farkas' lemma: :math:`Q\subseteq P` if and only if the following LP is feasible:
.. math::
\begin{align*}
\min_{\lambda, \mu, \mu_e} &\quad 0 \\
\text{s.t.} &\quad B = \lambda A + \mu A_e \\
&\quad B_e = \mu_e A_e
\end{align*}
"""
if isinstance(inner, np.ndarray):
# Testing if a single point is contained in the polyhedron.
if inner.ndim != 1 or inner.shape[0] != self.dim:
raise mpterr.DimMismatchError(self.dim, inner.shape[0])
if self.is_empty:
return False
if self.hasHrep:
return (np.all(self.A @ inner <= self.b + tolerances['zero']) and
np.all(np.abs(self.Ae @ inner - self.be) <= tolerances['zero']))
return self.distance(inner) <= tolerances['zero']
elif isinstance(inner, Polyhedron):
# Testing if one polyhedron is contained in another.
if inner.dim != self.dim:
raise mpterr.DimMismatchError(self.dim, inner.dim)
# An empty set is always contained in any other set
if inner.is_empty:
return True
# An empty set cannot contain a non-empty set (note that at this
# point, due to the previous check, we know that inner is not empty)
if self.is_empty:
return False
# Case: Inner is in V-representation, outer is in H- or V-representation
# Test that the rays of the inner polyhedron are contained in the outer polyhedron
if inner.hasVrep:
outer = Polyhedron(self.Hrep.homogenize)
homo = inner.Vrep.homogenize
return all(r in outer for r in homo.R)
# Inner is in H-rep - use Farkas' lemma.
# We need the outer to have an H-rep.
self._getHrep()
# Apply Farkas lemma to the homogenization
# inner := {x | A * x <= 0, Ae * x == 0}
# outer := {x | B * x <= 0, Be * x == 0}
#
# inner subseteq outer if and only if every inequality of outer can be explained by
# a (non-negative) combination of inner constraints
# B = Lambda @ A + Mu @ Ae
# Be = Mue @ Ae
# Lambda >= 0
i = inner.Hrep.homogenize
o = self.Hrep.homogenize
# Lambda = cp.Variable((o.A.shape[0], i.A.shape[0]))
# ineq = o.A - Lambda @ i.A
# if i.Ae.shape[0] > 0:
# Mu_ineq = cp.Variable((o.A.shape[0], i.Ae.shape[0]))
# ineq -= Mu_ineq @ i.Ae
# obj = cp.sum_squares(ineq)
#
# if o.Ae.shape[0] > 0:
# if i.Ae.shape[0] == 0:
# return False # Can't have an affine set in the outer without also an affine set in the inner
# Mu_eq = cp.Variable((o.Ae.shape[0], i.Ae.shape[0]))
# eq = o.Ae - Mu_eq @ i.Ae
# obj += cp.sum_squares(eq)
#
# con: List[cp.Constraint] = [Lambda[:] >= 0]
# prob = cp.Problem(cp.Minimize(obj), con)
# prob.solve(cp.GUROBI)
# if prob.status == cp.OPTIMAL:
# assert isinstance(prob.value, float), "cvxpy objective isn't a float?!"
# return float(prob.value) < tolerances['zero']
# elif prob.status == cp.INFEASIBLE:
# raise ValueError("This problem should never be infeasible... likely a CVXPY bug")
# elif prob.status == cp.UNBOUNDED:
# raise ValueError("This problem should never be unbounded... likely a CVXPY bug")
# else:
# raise RuntimeError
Lambda = cp.Variable((o.A.shape[0], i.A.shape[0]), nonneg=True)
constraints = []
# outer inequalities must be representable
if i.Ae.shape[0] > 0:
Mu_ineq = cp.Variable((o.A.shape[0], i.Ae.shape[0]))
constraints.append(Lambda @ i.A + Mu_ineq @ i.Ae - o.A <= tolerances['zero'])
constraints.append(Lambda @ i.A + Mu_ineq @ i.Ae - o.A >= -tolerances['zero'])
else:
constraints.append(Lambda @ i.A - o.A <= 1E-8)
constraints.append(Lambda @ i.A - o.A >= -1E-8)
# outer equalities handling
if o.Ae.shape[0] > 0:
if i.Ae.shape[0] == 0:
return False
Mu_eq = cp.Variable((o.Ae.shape[0], i.Ae.shape[0]))
constraints.append(Mu_eq @ i.Ae - o.Ae <= tolerances['zero'])
constraints.append(Mu_eq @ i.Ae - o.Ae >= -tolerances['zero'])
# feasibility problem
prob = cp.Problem(cp.Minimize(cp.Constant(0.)), constraints)
prob.solve() # active-set solvers like DAQP seems to perform better?
if prob.status in [cp.OPTIMAL, cp.OPTIMAL_INACCURATE]:
return True
elif prob.status in [cp.INFEASIBLE, cp.INFEASIBLE_INACCURATE]:
return False
else:
# should not happen
raise ValueError("CVXPY returned status " + prob.status)
elif isinstance(inner, ConvexSet):
# Check if the given non-polyhedral set is contained in the current set
if inner.dim != self.dim:
raise mpterr.DimMismatchError(self.dim, inner.dim)
# Test that inner is in our AffineHull
for a, b in zip(self.Ae, self.be):
if inner.support(a) - b > tolerances['zero']:
return False
if inner.support(-a) + b > tolerances['zero']:
return False
# Test if the support of each inequality of self is non-negative wrt the given set
for a, b in zip(self.A, self.b):
if inner.support(a) - b > tolerances['zero']:
return False
return True
else:
raise TypeError("Invalid type. Supported types are Vector, Polyhedron, and ConvexSet.")
def __eq__(self, other: object) -> bool:
"""
Return True if the two polyhedra are equivalent and False otherwise.
"""
if not isinstance(other, Polyhedron):
return False
if self.dim != other.dim:
return False
return self in other and other in self
def __neg__(self) -> "Polyhedron":
r"""
Unitary minus for a polyhedron.
For a Polyhedron with H-representation :math:`P = \{ x\in\mathbb{R}^n \mid Ax\leq b, A_ex=b_e \}`,
its unitary minus :math:`Q=-P` is defined as:
.. math::
Q = \{ x\in\mathbb{R}^n \mid -A x \leq b, -A_e x = b_e \}
For a Polyhedron with V-representation
:math:`P = \{ x\in\mathbb{R}^n \mid x = V^\top v + R^\top r, v,r \geq 0, 1^\top v = 1 \}`,
its unitary minus :math:`Q=-P` is defined as:
.. math::
Q = \{ x\in\mathbb{R}^n \mid x = -V^\top v - R^\top r, v,r \geq 0, 1^\top v = 1 \}
"""
if self.hasHrep and not self.hasVrep:
H_new = HData(A=-self.A, b=self.b, Ae=-self.Ae, be=self.be, is_minimal=self.Hrep.is_minimal)
return Polyhedron(H=H_new)
elif self.hasVrep and not self.hasHrep:
V_new = VData(V=-self.V, R=-self.R, is_minimal=self.Vrep.is_minimal)
return Polyhedron(V=V_new)
else:
H_new = HData(A=-self.A, b=self.b, Ae=-self.Ae, be=self.be, is_minimal=self.Hrep.is_minimal)
V_new = VData(V=-self.V, R=-self.R, is_minimal=self.Vrep.is_minimal)
return Polyhedron(H=H_new, V=V_new)
def __add__(self, other: Union[Vector, "Polyhedron"]) -> "Polyhedron":
"""
Compute the Minkowski sum of the polyhedron with a point or another polyhedron.
"""
return self.minkowski_sum(other)
[docs]
def minkowski_sum(self, other: Union[Vector, "Polyhedron"]) -> "Polyhedron":
r"""
Compute the Minkowski sum of the polyhedron with a point or another polyhedron.
- The Minkowski sum of a polyhedron :math:`P\subseteq\mathbb{R}^n` and a point :math:`z\in\mathbb{R}^n` is defined as:
.. math::
P \oplus z = \{ x+z \mid x\in P \}.
- The Minkowski sum of two polyhedra :math:`P\subseteq\mathbb{R}^n` and :math:`Q\subseteq\mathbb{R}^n` is defined as:
.. math::
P \oplus Q = \{ x+y \mid x\in P, y\in Q \}.
"""
if isinstance(other, np.ndarray):
if other.ndim != 1:
raise ValueError("The point must be a vector.")
if other.size != self.dim:
raise mpterr.DimMismatchError(self.dim, other.size)
# TODO: any better idea than performing vertex enumeration?
if not self.hasVrep:
self._getVrep()
vertices = self.V
# if the V-rep has no vertices, we should add a vertex (0,...,0) to sum with the point
if vertices.shape[0] == 0:
vertices = np.zeros((1, self.dim))
vertices_sum = vertices + np.tile(other, (vertices.shape[0], 1))
return Polyhedron(V=VData(V=vertices_sum, R=self.R))
elif isinstance(other, Polyhedron):
if other.dim != self.dim:
raise mpterr.DimMismatchError(self.dim, other.dim)
# if both have V-representations, just take sums of vertices and rays
if self.hasVrep and other.hasVrep:
vertices_sum = np.array([vp + vq for vp in self.V for vq in other.V])
rays_sum = np.array([rp + rq for rp in self.R for rq in other.R])
if vertices_sum.size == 0:
vertices_sum = np.zeros((0, self.dim))
if rays_sum.size == 0:
rays_sum = np.zeros((0, self.dim))
return Polyhedron(V=VData(vertices_sum, rays_sum))
# if both have H-representations, compute projection (w, x, y) onto w
elif self.hasHrep and other.hasHrep:
self.minHrep(inplace=True)
other.minHrep(inplace=True)
A1 = np.column_stack((np.zeros_like(self.A), self.A, np.zeros_like(self.A))) # [0, A1, 0] [w; x; y] <= b1
A2 = np.column_stack((np.zeros_like(other.A), np.zeros_like(other.A), other.A)) # [0, 0, A2] [w; x; y] <= b2
A_new = np.vstack((A1, A2))
b_new = np.concatenate((self.b, other.b))
Ae1 = np.column_stack((np.zeros_like(self.Ae), self.Ae, np.zeros_like(self.Ae))) # [0, Ae1, 0] [w; x; y] = be1
Ae2 = np.column_stack((np.zeros_like(other.Ae), np.zeros_like(other.Ae), other.Ae)) # [0, 0, Ae2] [w; x; y] = be2
Ae3 = np.column_stack((-np.eye(self.dim), np.eye(self.dim), np.eye(self.dim))) # [-I, I, I] [w; x; y] = 0
Ae_new = np.vstack((Ae1, Ae2, Ae3))
be_new = np.concatenate((self.be, other.be, np.zeros((self.dim,))))
return Polyhedron(HData(A_new, b_new, Ae_new, be_new)).projection([*range(self.dim)])
# If one of them has H-rep and another has V-rep, compute the V-rep
else:
self._getVrep()
other._getVrep()
return self.minkowski_sum(other)
else:
raise TypeError("Can only compute Minkowski sum with another point or polyhedron.")
def __sub__(self, other: Union[Vector, "Polyhedron"]):
"""
Compute the Pontryagin difference between a polyhedron with a point or another polyhedron.
"""
return self.pontryagin_difference(other)
[docs]
def pontryagin_difference(self, other: Union[Vector, "Polyhedron"]) -> "Polyhedron":
r"""
Compute the Pontryagin difference between a polyhedron with a point or another polyhedron.
- The Pontryagin difference of a polyhedron :math:`P\subseteq\mathbb{R}^n` with a point :math:`z\in\mathbb{R}^n` is defined as:
.. math::
P \ominus z = \{ x-z \mid x\in P \} = P \oplus (-z).
- The Pontryagin difference of a polyhedron :math:`P\subseteq\mathbb{R}^n` with another polyhedron :math:`Q\subseteq\mathbb{R}^n` is defined as:
.. math::
P \ominus Q = \{ x\in P \mid \forall y\in Q, x+y \in P \}.
"""
if isinstance(other, np.ndarray):
return self.minkowski_sum(-other)
if isinstance(other, Polyhedron):
if self.dim != other.dim:
raise mpterr.DimMismatchError(self.dim, other.dim)
self.minHrep()
other.minHrep()
# The affine hull of the other polyhedron must be a subset of the affine hull of the self polyhedron,
# otherwise the pontryagin difference is an empty set
empty_set = Polyhedron(HData(A=np.zeros((1,self.dim)), b=np.array([-1])))
if np.vstack((self.He, other.He)).size and other.He.size:
if np.linalg.matrix_rank(np.vstack((self.He, other.He))) > np.linalg.matrix_rank(other.He):
return empty_set
# If the two polyhedra are equal, return empty set
if self == other:
return empty_set
# Subtract the support of other for each inequality of P
b_new = np.zeros_like(self.b)
for i in range(self.A.shape[0]):
b_new[i] = self.b[i] - other.support(self.A[i, :])
if np.any(b_new == -np.inf):
return empty_set
# Remove remaining +Inf rows in b
inf_rows = (b_new == np.inf)
Hrep = HData(A=self.A[~inf_rows], b=b_new[~inf_rows], Ae=self.Ae, be=self.be)
return Polyhedron(Hrep)
[docs]
def cartesian_product(self, other: "Polyhedron") -> "Polyhedron":
r"""
Compute the Cartesian product with another polyhedron. The Cartesian product of two polyhedra
:math:`P\subseteq\mathbb{R}^m` and :math:`Q\subseteq\mathbb{R}^n` is defined as:
.. math::
P \times Q := \{ (x,y) \in\mathbb{R}^{m+n} \mid x\in P, y \in Q \}.
"""
if self.hasHrep and other.hasHrep:
A_new = np.vstack((
np.hstack((self.A, np.zeros((self.A.shape[0], other.dim)))),
np.hstack((np.zeros((other.A.shape[0], self.dim)), other.A))
))
b_new = np.concatenate((self.b, other.b))
Ae_new = np.vstack((
np.hstack((self.Ae, np.zeros((self.Ae.shape[0], other.dim)))),
np.hstack((np.zeros((other.Ae.shape[0], self.dim)), other.Ae))
))
be_new = np.concatenate((self.be, other.be))
return Polyhedron(H=HData(A=A_new, b=b_new, Ae=Ae_new, be=be_new))
elif self.hasVrep and other.hasVrep:
# Compute vertex pairs using broadcasting
vertices_product = np.hstack((
np.repeat(self.V, other.V.shape[0], axis=0),
np.tile(other.V, (self.V.shape[0], 1))
))
# Compute rays: stack rays from both polyhedra with zero-padding
rays_product = np.vstack((
np.hstack((self.R, np.zeros((self.R.shape[0], other.R.shape[1])))),
np.hstack((np.zeros((other.R.shape[0], self.V.shape[1])), other.R))
))
vrep = VData(V=vertices_product, R=rays_product)
return Polyhedron(V=vrep)
else:
raise NotImplementedError
def __mul__(self, other: Union[Scalar, "Polyhedron"]) -> "Polyhedron":
if isinstance(other, Polyhedron):
return self.cartesian_product(other)
elif np.isscalar(other):
T = float(other) * np.eye(self.dim)
return self.affine_map(T)
else:
raise TypeError("Unsupported type. Supported types are scalar and Polyhedron.")
__rmul__ = __mul__
__array_priority__ = 1.0 # higher than ndarray’s (~0.0), make A @ P work when A is a NumPy array
def __rmatmul__(self, T: Matrix) -> "Polyhedron":
if (not isinstance(T, np.ndarray)) or np.ndim(T) != 2:
raise TypeError("The transformation matrix T must be a matrix.")
return self.affine_map(T)
[docs]
def affine_map(self, T: Matrix, t: Optional[Vector] = None) -> "Polyhedron":
r"""
Computes an affine map :math:`P \rightarrow Q` of polyhedron :math:`P` to polyhedron :math:`Q`
based on the transformation matrix :math:`T`. The polyhedron :math:`Q` is given by
.. math::
Q = \{ y \in \mathbb{R}^n \mid y=Tx+t, x \in P \subseteq \mathbb{R}^d \}
The matrix :math:`T` must be real with :math:`n` rows and :math:`d` columns.
- If :math:`n < d` then this operation is referred to as projection.
- If :math:`n = d` then this operation is referred to as rotation/skew.
- If :math:`n > d` then this operation is referred to as lifting.
.. note::
This operation is also overloaded via ``__matmul__`` and ``__add__``. So for a Polyhedron object ``P``, ``P.affine_map(T, t)`` is equivalent to ``T @ P + t``.
.. seealso::
- :func:`inv_affine_map` - Compute the inverse affine map
"""
if np.ndim(T) != 2:
raise ValueError("The transformation matrix T must be a matrix.")
if not T.shape[1] == self.dim:
raise ValueError("The transformation matrix T must be have {} columns.".format(self.dim))
t = np.zeros((T.shape[0],)) if t is None else np.reshape(t, (-1,))
if t.size != T.shape[0]:
raise ValueError("The transformation matrix T and translation vector t must have same number of rows.")
n, d = T.shape[0], T.shape[1]
# The affine map of an empty set is still an empty set
if self.is_empty:
return Polyhedron.empty_set(n)
# If T is (close to) zero matrix, the image is a single point t
if np.linalg.matrix_norm(T, ord=2) < tolerances['zero']:
return Polyhedron.singleton(t)
Hrep_new, Vrep_new = None, None
if self.hasVrep:
T_transpose = T.transpose()
t_transpose = np.reshape(t, (1, -1))
V_new = self.V @ T_transpose + t_transpose
R_new = self.R @ T_transpose
Vrep_new = VData(V=V_new, R=R_new)
if self.hasHrep:
self._compute_affine_dimension_for_Hrep() # extract implicit equalities
## T is invertible, leads to a full-dimensional image
if n == d and abs(np.linalg.det(T)) > tolerances['zero']:
# x = T^-1 * (y-t) => A x<=b, Ae x = be
T_inv = np.linalg.inv(T)
Hrep_new = HData(A=self.A @ T_inv, b=self.b + self.A @ T_inv @ t, Ae=self.Ae @ T_inv, be=self.be + self.Ae @ T_inv @ t)
# T is singular or rectangular, leading to a lower-dimensional image
# Q = Proj_y {(x,y) | Ax <= b, Ae x = be, y = Tx + t}
else:
# [A, -I] [x; y] <= b
A_xy = np.hstack((self.A, -np.zeros((self.A.shape[0], n))))
b_xy = self.b
Ae1 = np.hstack((self.Ae, np.zeros((self.Ae.shape[0], n)))) # [Ae, 0] [x; y] = be
Ae2 = np.hstack((-T, np.eye(n))) # [-T, I] [x; y] = t
Ae_xy = np.vstack((Ae1, Ae2))
be_xy = np.concatenate((self.be, t))
Q_xy = Polyhedron.from_Hrep(A=A_xy, b=b_xy, Ae=Ae_xy, be=be_xy)
Q = Q_xy.projection(list(range(d, n + d)))
Hrep_new = Q.Hrep
return Polyhedron(H=Hrep_new, V=Vrep_new)
[docs]
def inv_affine_map(self, T: Matrix, t: Optional[Vector] = None) -> "Polyhedron":
r"""
Computes an inverse affine map :math:`Q` of polyhedron :math:`P` based on the transformation matrix :math:`T`
and vector :math:`t`. The polyhedron :math:`Q` is given by:
.. math::
Q = \{ x \in \mathbb{R}^n \mid T x + t \in P \subseteq \mathbb{R}^n \}
The matrix :math:`T` must be a square real matrix. The vector :math:`t`, if omitted, defaults to a zero vector
of corresponding dimension.
.. seealso::
- :func:`affine_map` - Compute the affine map
"""
if t is None:
t = np.zeros((self.dim,))
raise NotImplementedError
@property
def affine_dimension(self) -> int:
r"""
Compute the affine dimension of the polyhedron, i.e.,
the number of degrees of freedom you have when moving within the affine set.
.. admonition:: Math
The implementation is based on the method proposed in Section 7.3 of
`Fukuda's Polyhedral Computation <https://www.research-collection.ethz.ch/handle/20.500.11850/426218>`_.
- For a V-polyhedron :math:`P = \{ x \in \mathbb{R}^n \mid x = V^\top \lambda + R^\top\mu, 1^\top \lambda = 1, \mu \geq 0, \lambda \geq 0\}`,
.. math::
\text{affdim}(P) = \text{rank} \begin{bmatrix} V^\top & R^\top \\ 1_{1\times n} & 0_{1\times n} \end{bmatrix} -1.
- For a H-polyhedron :math:`P = \{ x \in \mathbb{R}^n \mid A x \leq b, A_e x = b_e \}`:
.. math::
\text{affdim}(P) = n - \text{rank}\tilde{A}_e.
where :math:`\tilde{A}_e` is the matrix containing all explicit and implicit equalities.
.. Note::
If a polyhedron is empty, its affine dimension is defined as -1.
"""
# If already computed, directly return the result
if self._affine_dim is np.nan:
if self.hasVrep:
Vrep_homo = self._Vrep.homogenize
self._affine_dim = np.linalg.matrix_rank(Vrep_homo.R, tol=tolerances['zero']) - 1
elif self.hasHrep:
self._compute_affine_dimension_for_Hrep()
else:
raise ValueError("Neither H-representation nor V-representation found")
return self._affine_dim
def _compute_affine_dimension_for_Hrep(self) -> None:
r"""
Determine if the H-polyhedron is full-dimensional, empty, or lower-dimensional by solving one single LP.
The implementation is based on (7.4) in
`Fukuda's Polyhedral Computation <https://www.research-collection.ethz.ch/handle/20.500.11850/426218>`_.
This function modifies self._affine_dim if the polyhedron is either full-dim or empty.
If the polyhedron is lower-dimensional, self._affine_dim is still None
"""
if not self.hasHrep:
raise RuntimeError("This function is for H-representation only because we want to avoid doing facet enumeration.")
# Check if the equalities result in an empty set
if self.Ae.size != 0:
sol = np.linalg.lstsq(self.Ae, self.be)[0]
if not np.allclose(self.Ae @ sol, self.be, atol=tolerances['zero']):
self._affine_dim = -1
# Solve an LP to determine if it is full-dimensional
# This LP is always feasible and bounded
x0 = cp.Variable()
x = cp.Variable(self.dim)
constraints = []
if self.A.size:
constraints = [self.A @ x + np.ones(self.A.shape[0]) * x0 <= self.b]
if self.Ae.size:
constraints += [self.Ae @ x + np.ones(self.Ae.shape[0]) * x0 <= self.be]
constraints += [-self.Ae @ x + np.ones(self.Ae.shape[0]) * x0 <= -self.be]
constraints += [x0 <= 1]
prob = cp.Problem(cp.Maximize(x0), constraints)
prob.solve()
if prob.status != cp.OPTIMAL:
raise RuntimeError("The LP when determining if the problem is full dim is not solved to optimality.")
# full dimensional
if prob.value >= tolerances['zero']:
self._affine_dim = self.dim
# empty set
elif prob.value <= -tolerances['zero']:
self._affine_dim = -1
# lower dimensional
else:
# Detect implicit equalities
implicit_equalities = np.abs(self.A @ x.value - self.b) <= tolerances['zero']
Ae_new = np.vstack((self.Ae, self.A[implicit_equalities, :]))
be_new = np.concatenate([self.be, self.b[implicit_equalities]])
A_new = self.A[np.logical_not(implicit_equalities), :]
b_new = self.b[np.logical_not(implicit_equalities)]
self._affine_dim = self.dim - np.linalg.matrix_rank(Ae_new)
# TODO: Should we make the implicit equalities into explicit or not?
self._Hrep = HData(A=A_new, b=b_new, Ae=Ae_new, be=be_new)
@property
def codimension(self) -> int:
r"""
Returns the codimension of the polyhedron, i.e., the number of constraints/dimensions
by which the affine set is "lower" than the ambient space.
For example, suppose we have a non-empty polyhedron
.. math::
P = \{ x \in \mathbb{R}^n \mid A x \leq b, A_e x = b_e \},
where the rank of :math:`A_e` (number of rows) is :math:`m`.
The codimension of :math:`P` is :math:`m`.
"""
# TODO: should we define the codimension of an empty set as infinity?
return self.dim - self.affine_dimension
[docs]
def plot(self, ax: Union[Axes, Plotter], **kwargs):
r"""
Plot the Polyhedron.
.. Note::
Require V-representation. Will be computed if necessary.
"""
if self.dim > 3:
raise ValueError("Can only plot 2D or 3D polyhedra")
if self.dim <= 1:
raise ValueError("Cannot plot 1D or 0D polyhedra")
if isinstance(ax, Axes):
from mpt4py.geometry.visualization.plot_matplotlib import MatplotlibPlotter
plotter_protocal = MatplotlibPlotter(ax)
plotter_protocal.plot_convexhull(self.Vrep.V, self.Vrep.R, **kwargs)
elif isinstance(ax, Plotter):
from mpt4py.geometry.visualization.plot_pyvista import PyvistaPlotter
plotter_protocal = PyvistaPlotter(ax)
plotter_protocal.plot_convexhull(self.Vrep.V, self.Vrep.R, **kwargs)
else:
raise NotImplementedError("Unsupported backend. Supported backends are matplotlib and pyvista.")
[docs]
def fplot(self, ax: Union[Axes, Plotter], func_name: Optional[str] = None, **kwargs):
r"""
Plot a single function over a polyhedron.
.. Note::
Require H-representation. Will be computed if necessary.
"""
if self.dim != 2:
raise ValueError("Can only do fplot for 2D polyhedra")
# If no function name is given, and there is only one function, use that one
if len(self._functions) == 1 and func_name is None:
func_name = list(self._functions.keys())[0]
func = self.get_function(func_name)
from mpt4py.functions.affine_function import AffineFunction
if isinstance(func, AffineFunction):
vertices_xy = self.V
vertices_z = np.zeros((vertices_xy.shape[0], 1))
rays_xy = self.R
rays_z = np.zeros((rays_xy.shape[0], 1))
for i in range(vertices_xy.shape[0]):
vertices_z[i] = func.evaluate(vertices_xy[i, :])
for i in range(rays_xy.shape[0]):
rays_z[i] = func.evaluate(rays_xy[i, :])
vertices = np.hstack((vertices_xy, vertices_z.reshape(-1, 1)))
rays = np.hstack((rays_xy, rays_z.reshape(-1, 1)))
if isinstance(ax, Axes):
from mpt4py.geometry.visualization.plot_matplotlib import MatplotlibPlotter
plotter_protocal = MatplotlibPlotter(ax)
elif isinstance(ax, Plotter):
from mpt4py.geometry.visualization.plot_pyvista import PyvistaPlotter
plotter_protocal = PyvistaPlotter(ax)
else:
raise NotImplementedError("Unsupported backend. Supported backends are matplotlib and pyvista.")
plotter_protocal.plot_convexhull(points=vertices, rays=rays, **kwargs)
else:
raise NotImplementedError("Only affine functions are supported for fplot currently.")
[docs]
def lift(self, num_dims: int) -> "Polyhedron":
r"""
Lift the polyhedron :math:`P\in\mathbb{R}^n` to a higher dimension in :math:`\mathbb{R}^{n+m}`:
.. math::
\{ (x, y) \mid x \in P, y \in \mathbb{R}^m \}.
"""
if num_dims < 1:
raise ValueError("Invalid number of dimensions")
args = {}
if self.hasHrep:
Ho = self.Hrep
args['H'] = HData(A=np.hstack([Ho.A, np.zeros((Ho.A.shape[0], num_dims))]), b=Ho.b,
Ae=np.hstack([Ho.Ae, np.zeros((Ho.Ae.shape[0], num_dims))]), be=Ho.be)
if self.hasVrep:
Vo = self.Vrep
args['V'] = VData(V=np.hstack([Vo.V, np.zeros((Vo.V.shape[0], num_dims))]),
R=np.hstack([Vo.R, np.zeros((Vo.R.shape[0], num_dims))]))
return Polyhedron(**args)
[docs]
def intersect(self, other: "Polyhedron") -> "Polyhedron":
"""
Compute the intersection of this polyhedron with another Polyhedron.
"""
# Requires the conversion to H-representation
if self.dim != other.dim:
raise mpterr.GeometryError('Cannot intersect polyhedra of different dimensions')
H = HData(A=np.vstack([self.A, other.A]), b=np.hstack([self.b, other.b]),
Ae=np.vstack([self.Ae, other.Ae]), be=np.hstack([self.be, other.be]))
return Polyhedron(H)
[docs]
def homogenize(self) -> "Polyhedron":
r"""
Homogenize the Polyhedron into a cone.
- For an H-representation :math:`P = \{ x\in\mathbb{R}^n \mid A x \leq b, A_e x = b_e \}`, its homogenization is:
.. math::
\bar{P} = \{ y\in\mathbb{R}^{n+1} \mid [A ~-b] \leq 0, [A_e ~-b_e] x = 0 \}.
- For a V-representation: :math:`P = \{ x\in\mathbb{R}^n \mid x = V^\top v + R^\top r, r \geq 0, v\geq 0, 1^\top v = 1 \}`, its homogenization is:
.. math::
P_h = \{ y\in\mathbb{R}^{n+1} \mid y=\begin{bmatrix} V^\top & R^\top \\ 1 & 0\end{bmatrix} l, l \geq 0 \}.
"""
args = {}
if self.hasHrep:
args['H'] = self.Hrep.homogenize
if self.hasVrep:
args['V'] = self.Vrep.homogenize
return Polyhedron(**args)
@property
def is_bounded(self):
r"""
Returns `True` if the polyhedron is bounded and `False` otherwise.
"""
if self.hasVrep:
return self.R.shape[0] == 0
# We're unbounded if we can find a vector d such that A @ d <= 0, Ae @ d = 0
# Why? Let x be in the Polyhedron. Then x + d * alpha is also in the polyhedron for any alpha >= 0.
# Assume A @ x <= b, Ae @ x = be. Then A @ (x + d * alpha) = A @ x + A @ d * alpha <= b
var = cp.Variable(self.dim)
c = cp.Parameter(self.dim)
con: List[cp.Constraint] = [self.A @ var <= 0]
con += [-1 <= var, var <= 1] # Add box constraints to avoid numerical issues with unbounded polyhedra
if self.Ae.size > 0:
con += [self.Ae @ var == 0]
prob = cp.Problem(cp.Maximize(c.T @ var), con)
for i in range(self.dim):
c.value = np.zeros(self.dim)
c.value[i] = 1
prob.solve()
if prob.status != cp.OPTIMAL:
raise ValueError("Error solving LP to find an unbounded direction")
if float(prob.value) > tolerances['zero']:
return False
c.value = np.zeros(self.dim)
c.value[i] = -1
prob.solve()
if prob.status != cp.OPTIMAL:
raise ValueError("Error solving LP to find an unbounded direction")
if float(prob.value) > tolerances['zero']:
return False
return True
@property
def is_full_dim(self) -> bool:
r"""
Determine if the polyhedron is full dimensional or not.
A polyhedron :math:`P\in\mathbb{R}^n` is full dimensional if and only if the dimension of its affine dimension
is the same as :math:`n` (its dimension of representation), or equivalently if there exists a non-empty ball
of dimension :math:`n` that is contained within it.
"""
return self.affine_dimension == self.dim
@property
def is_empty(self) -> bool:
r"""
Returns `True` if the polyhedron is empty and `False` otherwise.
"""
return self.affine_dimension == -1
[docs]
def affine_hull(self, inplace: bool = False) -> AffineSet:
r"""
Compute the affine hull of the polyhedron, which is the smallest affine set
that contains it.
This function detects if the polyhedron has an implicit affine set,
e.g., if the inequalities contain :math:`\alpha^\top x \leq \beta` and :math:`\alpha^\top x \geq \beta`.
If `inplace == True`, then the HData of this object is modified.
"""
if self.hasVrep:
homo = self.homogenize()
# Affine hull is x = homo.R' * r
aff = AffineSet.from_basis(homo.R.T, homo.V[0])
# Un-lift the result A[x;-1] = 0
return AffineSet(aff.A[:,:-1], aff.A[:,-1])
if self.hasHrep:
# For each inequality, we want to know if it's part of the affine hull.
#
# Test if Ai x = c is in the affine hull
#
# Ji = min Ai s.t. x in P
#
# Ji == c if and only if Ai x = c is in the affine hull
is_affine = [abs(-self.support(-a) - b) < tolerances['zero'] for a, b in zip(self.A, self.b)]
Ae = np.vstack([self.Ae, self.A[is_affine]])
be = np.hstack([self.be, self.b[is_affine]])
aff = AffineSet(A=Ae, b=be).minimal_representation()
if inplace:
assert isinstance(self._Hrep, HData), "Invalid H-representation" # for pylance
self._Hrep.Ae = aff.A
self._Hrep.be = aff.b
self._Hrep.A = self.A[not is_affine]
self._Hrep.b = self.b[not is_affine]
return aff
raise ValueError("Invalid representation - should never hit this point")
[docs]
@staticmethod
def bounding_box(obj: ConvexSet) -> "Polyhedron":
r"""
Compute the minimum bounding box of the polyhedron.
If the polyhedron is unbounded, then the box will be unbounded in those directions.
Note that if the polyhedron has an affine hull, then the bounding box will
take the form
.. math::
B = \{ x\in\mathbb{R}^n \mid x_{lb} \leq x \leq x_{ub}, A_e x = b_e \}.
"""
A = []
b = []
for i in range(obj.dim):
c = np.zeros(obj.dim)
c[i] = 1
sup = obj.support(c)
if sup is not np.inf:
A.append(c)
b.append(sup)
sup = obj.support(-c)
if sup is not np.inf:
A.append(-c)
b.append(sup)
A = np.array(A)
b = np.array(b)
try:
Ae, be = obj.affine_hull().coefficients
return Polyhedron(HData(A=A, b=b, Ae=Ae, be=be))
except AttributeError:
return Polyhedron(HData(A=A, b=b))
[docs]
def cheby_center(self, raise_on_empty_set: Optional[bool] = True) -> Tuple[Vector, float]:
r"""
Compute the center and radius of the polyhedron's Chebyshev ball.
Given a polyhedron with H-representation, the center :math:`c` and radius :math:`r` of the Chebyshev ball
is computed by solving the following Linear Programming:
.. math::
\begin{align*}
\max_{c,r} &\quad r \\
\text{s.t.} &\quad a_i c + \|a_i\|_2 r \leq b_i, \quad \forall i\in \{1,...,m\} \\
&\quad A_e c = b_e
\end{align*}
where :math:`a_i` is the :math:`i` th row of :math:`A` and :math:`m` is the number of rows of :math:`A`.
.. note::
- If the polyhedron is empty, the Chebyshev center is set to :math:`-\infty` and the radius is set to :math:`-\infty`.
- If the polyhedron is unbounded, the Chebyshev center is set to :math:`+\infty` and the radius is set to :math:`+\infty`.
.. warning::
- If the polyhedron is lower-dimensional and it only has a V-representation, then its H-representation is non-unique. Different H-representations give different Chebyshev balls.
"""
if np.any(np.isnan(self._chebyshev_center)) or self._chebyshev_radius is np.nan:
self._cheby_center()
if self._chebyshev_radius > tolerances['zero']:
return self._chebyshev_center, self._chebyshev_radius
elif self._chebyshev_radius > -tolerances['zero']:
return self._chebyshev_center, 0.
else:
if raise_on_empty_set:
raise mpterr.InfeasibleError("The polyhedron is empty!")
else:
return self._chebyshev_center, self._chebyshev_radius
def _cheby_center(self) -> None:
"""
Compute the Chebyshev center.
"""
if not self.hasHrep:
self._getHrep()
c = cp.Variable(self.dim) # center of the chebyshev ball
r = cp.Variable() # radius of the chebyshev ball
constraints = []
if self.b.size:
constraints += [self.A @ c + np.linalg.norm(self.A, axis=1) * r <= self.b]
if self.be.size:
constraints += [self.Ae @ c == self.be]
objective = cp.Maximize(r)
problem = cp.Problem(objective, constraints)
problem.solve(solver=cp.ECOS, abstol=tolerances['zero'])
if problem.status == cp.UNBOUNDED:
chebyshev_center = np.inf * np.ones(self.dim)
chebyshev_radius = np.inf
elif problem.status == cp.OPTIMAL:
chebyshev_center = c.value # type: ignore
chebyshev_radius = r.value # type: ignore
elif problem.status == cp.INFEASIBLE:
chebyshev_center = -np.inf * np.ones(self.dim)
chebyshev_radius = -np.inf
else:
raise RuntimeError("Error in solving LP")
# obj_quad_coeff = np.zeros((self.dim + 1, self.dim + 1))
# obj_lin_coeff = np.zeros(self.dim + 1)
# obj_lin_coeff[-1] = -1.
# ineq_lhs = np.hstack((self.A, np.linalg.norm(self.A, axis=1).reshape(-1, 1))) if self.b.size else None
# ineq_rhs = self.b if self.b.size else None
# eq_lhs = np.hstack((self.Ae, np.zeros((self.Ae.shape[0], 1)))) if self.be.size else None
# eq_rhs = self.be if self.be.size else None
#
# solver = piqp.DenseSolver()
# solver.settings.verbose = False
# solver.settings.compute_timings = False
# solver.settings.max_iter = 500
# solver.setup(P=obj_quad_coeff, c=obj_lin_coeff, A=eq_lhs, b=eq_rhs, G=ineq_lhs, h=ineq_rhs)
# status = solver.solve()
# if status == piqp.Status.PIQP_DUAL_INFEASIBLE:
# chebyshev_center = np.inf * np.ones(self.dim)
# chebyshev_radius = np.inf
# elif status == piqp.Status.PIQP_SOLVED:
# chebyshev_center = solver.result.x[:-1]
# chebyshev_radius = solver.result.x[-1]
# elif status == piqp.Status.PIQP_PRIMAL_INFEASIBLE:
# chebyshev_center = -np.inf * np.ones(self.dim)
# chebyshev_radius = -np.inf
# else:
# raise RuntimeError("Error in solving LP")
# # Call scipy.optimize.linprog
# sol = linprog(c=obj_coeff, A_ub=ineq_lhs, b_ub=ineq_rhs, A_eq=eq_lhs, b_eq=eq_rhs, method='highs')
# if sol.status == 0:
# chebyshev_center = sol.x[:-1]
# chebyshev_radius = sol.x[-1]
# elif sol.status == 2:
# chebyshev_center = -np.inf * np.ones(self.dim)
# chebyshev_radius = -np.inf
# elif sol.status == 3:
# chebyshev_center = np.inf * np.ones(self.dim)
# chebyshev_radius = np.inf
# else:
# raise RuntimeError("Error in solving LP")
self._chebyshev_center = chebyshev_center
self._chebyshev_radius = chebyshev_radius
# If the radius is strictly negative, this implies the polyhedron is an empty set
if self._chebyshev_radius <= tolerances["zero"]:
self._affine_dim = -1
[docs]
def interior_point(self, facet_index: Optional[int] = None) -> Tuple[Vector, bool]:
r"""
Compute an interior point of the polyhedron.
Args:
facet_index (int): the index of the facet to compute the interior point for.
If None, compute the interior point of the whole polyhedron.
Returns:
An interior point of the polyhedron and a boolean indicating whether the point is strictly inside the polyhedron.
"""
if self.is_empty:
raise mpterr.InfeasibleError("The polyhedron is empty.")
if facet_index is None:
if self.hasVrep:
point = np.mean(self.V, axis=0) if self.V.size > 0 else np.zeros(self.dim)
if self.R.size:
normalized_rays = self.R / np.linalg.norm(self.R, axis=1, keepdims=True)
point = np.mean(normalized_rays, axis=0)
# determine if the point is strictly inside the polyhedron
strict_interior = self.is_full_dim
return point, strict_interior
else:
# if the polyhedron is an affine set, just compute x = Ae \ be
if self.A.size == 0:
return np.linalg.solve(self.Ae, self.be), False
x = cp.Variable(self.dim)
constraints = []
constraints += [self.A @ x <= self.b] if self.b.size else []
constraints += [self.Ae @ x == self.be] if self.be.size else []
objective = cp.Maximize(0)
problem = cp.Problem(objective, constraints)
problem.solve()
if problem.status == cp.OPTIMAL:
# verify the result
strict_interior = np.all(self.A @ x.value <= self.b - tolerances['zero']) and self.Ae.size == 0
return x.value, strict_interior
elif problem.status == cp.INFEASIBLE:
self._affine_dim = -1
raise mpterr.InfeasibleError("The polyhedron is empty!")
else:
raise RuntimeError("Error in solving LP.")
else:
if not self.hasHrep or not self.Hrep.is_minimal:
raise RuntimeError("Polyhedron must be in minimal H-representation when you want compute any interior point of its facets. Call minHRep() first.")
if self.A.size == 0:
raise ValueError("The polyhedron is an affine set. It does not have facets.")
if facet_index not in range(self.A.shape[0]):
raise ValueError("Invalid facet index. A valid facet index is in the range [0, {}].".format(self.A.shape[0] - 1))
raise NotImplementedError
[docs]
def sample(self, n_samples: int, strictly_interior: bool = False,
method: Literal['hit_and_run', 'reject'] = 'reject',
*, rng: Optional[np.random.Generator] = None, batch_size: int = 1024,
max_attempts: Optional[int] = None, whitening: bool = True, verbose: bool = False) -> np.ndarray:
"""
Uniformuly sample random points inside the polyhedron. Requires the polyhedron to be bounded and non-empty.
"""
if not self.is_bounded:
raise mpterr.UnboundedError("Sampling requires the polyhedron to be bounded.")
if self.is_empty:
raise mpterr.InfeasibleError("Sampling requires the polyhedron to be non-empty.")
# Set up random number generator for reproducibility
rng = np.random.default_rng(rng)
max_attempts = max(10_000, 100 * n_samples) if max_attempts is None else max_attempts
assert max_attempts >= n_samples, "max_attempts must be no smaller than n_samples"
def _rejection_sampling(P: "Polyhedron", n_samples: int) -> np.ndarray:
# We assume P is full-dimensional and has no equality constraints
assert P.is_full_dim
assert P.Ae.shape[0] == 0 and P.be.size == 0
# Axis-aligned bounding box from vertices
bbox = Polyhedron.bounding_box(P) # TODO: cache bounding_box() if it is called often?
lb = np.min(bbox.V, axis=0)
ub = np.max(bbox.V, axis=0)
samples = np.empty((0, P.dim))
attempts = 0
while len(samples) < n_samples and attempts < max_attempts:
m = min(batch_size, n_samples - len(samples))
X = rng.uniform(lb, ub, size=(m, P.dim))
mask = np.all(P.A @ X.T <= (P.b + tolerances['zero'] * (-1 if strictly_interior else 1))[:, None], axis=0)
samples = np.vstack((samples, X[mask])) if np.any(mask) else samples
attempts += m
if len(samples) == 0:
raise RuntimeError(f"Rejection sampling failed: no accepted points after {max_attempts} attempts. Check bounding box or try hit_and_run.")
out = samples
if out.shape[0] < n_samples:
raise RuntimeError(
f"Rejection sampling produced {out.shape[0]} < {n_samples}. "
"Bounding box likely too loose; try hit_and_run."
)
if verbose:
print(f"Successfully sampled {n_samples} after {attempts} attempts.")
return out[:n_samples]
def _dikin_whitening(P: "Polyhedron") -> Tuple[Vector, Matrix]:
"""
Use Dikin ellipsoid to whiten the polyhedron, i.e., transform it into a nearly isotropic position.
By applying the affine transformation, a tighter bounding box can be computed, in which the accept rate in rejection sampling will be higher.
"""
# We assume P is full-dimensional and has no equality constraints
assert P.is_full_dim
assert P.Ae.shape[0] == 0 and P.be.size == 0
xc, _ = P.cheby_center()
s = P.b - P.A @ xc # slacks in the inequalities
if np.any(s <= tolerances['zero']):
raise RuntimeError("Interior point not strictly interior.")
W = 1.0 / s ** 2
H = (P.A.T * W) @ P.A # Dikin Hessian. if P is full-dim, H should be positive definite
return xc, H
if method == "reject":
# Detect all implicit equalities and put them into Ae*x=be
self._getHrep()
self._compute_affine_dimension_for_Hrep()
# deal with equality constraints
if self.affine_dimension < self.dim:
A_reduced, b_reduced = self.A, self.b
# find a special solution satisfying Ae*x0 = be
x0, residuals, rank, s = np.linalg.lstsq(self.Ae, self.be, rcond=None)
if np.linalg.norm(self.Ae @ x0 - self.be, np.inf) > tolerances['zero']:
raise mpterr.InfeasibleError("Equality constraints has no solution, indicating the Polyhedron is empty!")
# compute null space basis
N = null_space(self.Ae)
# if the polyhedron only contains one single point
if self.affine_dimension == 0:
if all(self.A @ x0 <= self.b):
return np.tile(x0, (n_samples, 1))
else:
raise mpterr.InfeasibleError("The polyhedron should be empty but not detected previously ...")
# transform ineualities into lower-dimensional space: A*(x0 + N*y) <= b
A_reduced = self.A @ N
b_reduced = self.b - self.A @ x0
# sample in the lower dimensional space (affine hull)
P_reduced = Polyhedron.from_Hrep(A_reduced, b_reduced)
assert P_reduced.is_bounded
P_reduced._affine_dim = self.affine_dimension # ! we know P_reduced is full-dim. So directly mark it to avoid computing the affine dimension again
samples_reduced = P_reduced.sample(n_samples, strictly_interior=strictly_interior, method='reject', rng=rng, batch_size=batch_size, max_attempts=max_attempts, whitening=whitening, verbose=verbose)
samples = x0 + samples_reduced @ N.T # map back to the original space: x = x0 + N*y
else:
# sample in the original space
if not whitening:
samples = _rejection_sampling(self, n_samples)
else:
# whiten the polyhedron
xc, dikin_hessian = _dikin_whitening(self)
_, eigvecs = np.linalg.eigh(dikin_hessian) # extract the main direction of the polyhedron
rotation_mat = eigvecs
# for all x s.t. A*x<=b, find y = R*(x-xc) => A*(R'*y+xc)<=b
P_whitened = Polyhedron.from_Hrep(self.A @ rotation_mat.T, self.b - self.A @ xc)
P_whitened._affine_dim = self.dim # ! we know it is full-dim. Directly mark it to avoid computing the affine dimension again
# sample in the transformed space
samples = _rejection_sampling(P_whitened, n_samples)
# map back from whitening: x = R'y + xc
samples = samples @ rotation_mat + xc
return samples
elif method == "hit_and_run":
raise NotImplementedError
else:
raise ValueError("Unknown sampling method: {}".format(method))
[docs]
def projection(self, dims: Sequence[int], method: Optional[str] = None) -> "Polyhedron":
"""
Compute the orthogonal projection of the polyhedron on given dimensions.
"""
if self.hasVrep:
return Polyhedron(V=VData(self.V[:, dims], self.R[:, dims]))
# TODO: implement auto-selecting computation methods
if method is None:
method = 'vrep'
if method == 'vrep':
Vrep = self.Vrep
return Polyhedron(V=VData(self.V[:, dims], self.R[:, dims]))
elif method == 'fourier':
# Notice the fourier elimination in cdd cannot handle linearity, i.e., cannot handle equality constraints
# hrep_new = cdd_block_elimination(self.Hrep, dims)
hrep_new = cdd_fourier_elimination(self.Hrep, dims)
return Polyhedron(H=hrep_new)
elif method == 'ifourier':
raise NotImplementedError
elif method == 'mplp':
raise NotImplementedError
else:
raise ValueError("Unsupported method! Supported methods are: 'vrep', 'fourier', 'ifourier', 'mplp'.")
@property
def volume(self) -> float:
r"""
Compute the volume of this polyhedron.
.. Note::
The volume of a polyhedron is defined as 0 if it is lower-dimensional or empty and :math:`\infty` if it is full-dimensional and unbounded.
"""
if not self.is_full_dim:
return 0.
elif not self.is_bounded:
return np.inf
else:
hull = ConvexHull(self.V)
return hull.volume
[docs]
def get_facet(self, index: Optional[Union[int, list[int]]] = None) -> Union["Polyhedron", list["Polyhedron"]]:
r"""
Extract the facets of the polyhedron specified by the inequality indices. Each facet
is a polyhedron formed by rewriting the indexed ineqality as an equality constraint.
For a polyhedron :math:`P` in H-representation:
.. math::
P = \{ x\in\mathbb{R}^n \mid a_i^\top x \leq b_i, i \in \{1, ..., m\},
A_e x = b_e \}
given an index :math:`k`, the facet polyhedron :math:`Q` is built as:
.. math::
Q = \{ x\in\mathbb{R}^n \mid a_i^\top x \leq b_i, i \in \{1, ..., m\} \backslash \{k\},
a_k^\top x = b_k, A_e x = b_e \}
Args:
index (int): an integer or list of integers. If not provided, a list of all facets will be returned.
Note:
The polyhedron must have minimal (irredundant) H-representation,
otherwise an error will be thrown. Consider calling :func:`minHrep` beforehand.
"""
# check if the polyhedron has a min H-rep.
if not self.hasHrep:
raise ValueError("Polyhedron must be in minimal H-representation to extract facets. Call `minHrep` first.")
if not self.Hrep.is_minimal:
raise ValueError("H-representation of polyhedron is not minimal. Call `minHrep` first.")
if index is None:
all_indices = [*range(self.b.size)]
return self.get_facet(all_indices)
if isinstance(index, int):
if index < 0 or index >= self.dim:
raise ValueError("Invalid index: should be between 0 and {}".format(self.dim))
A_new = np.delete(self.A, index, 0)
b_new = np.delete(self.b, index, 0)
Ae_new = np.vstack((self.A[index, :], self.Ae))
be_new = np.concatenate((self.b[[index]], self.be))
return Polyhedron(HData(A_new, b_new, Ae_new, be_new))
elif isinstance(index, list):
if not all(isinstance(i, int) for i in index):
raise ValueError("Invalid index type: should be int or list of int.")
facets = []
for i in index:
facets.append(self.get_facet(i))
return facets
else:
raise TypeError("Invalid index type: should be int or list of int.")
[docs]
def is_neighbor(self, other: "Polyhedron") -> Tuple[bool, Union[int, None], Union[int, None]]:
r"""Test if a polyhedron touches another polyhedron along a given facet.
Return `True` if the polyhedron :math:`P` shares with another polyhedron :math:`Q` a part of a common facet. Both
polyhedra must be in H-representation. If they are not, the irredundant H-representations will be computed.
The function tests if the intersection of two polyhedra :math:`P` and S in dimension d is nonempty
and is of dimension :math:`d-1`. If this holds, then the two polyhedra are touching themself along
one facet.
.. Note::
If one of the polyhedra is empty, the function will return `False`.
.. Seealso::
- :func:`is_adjacent` - Determine if two polyhedra are adjacent.
"""
if self.dim != other.dim:
raise ValueError("Polyhedra must have the same dimension to check for neighbors.")
if self.is_empty or other.is_empty:
return False, None, None
# Compute the intersection of the two polyhedra
intersection = self.intersect(other)
# Check if the intersection is non-empty and has dimension (d-1)
if intersection.affine_dimension != self.dim - 1:
return False, None, None
else:
raise NotImplementedError
[docs]
def is_adjacent(self, other: "Polyhedron") -> bool:
r"""Test if a polyhedron shares a facet with another polyhedron.
Return true if the polyhedron P has a facet-to-facet property with the polyhedron Q. Both
polyhedra must be in H-representation. If they are not, the irredundant H-representations
will be computed.
The function tests if polyhedra P and Q are adjacent by checking if their intersection is of
dimension :math:`d-1` and if the intersection is a facet of both polyhedra.
.. Seealso::
- :func:`is_neighbor` - Determine if two polyhedra are neighbors.
"""
if self.dim != other.dim:
raise ValueError("Polyhedra must have the same dimension to check for adjacency.")
# Ensure both polyhedra have minimal H-representations
self.minHrep(inplace=True)
other.minHrep(inplace=True)
# Compute the intersection of the two polyhedra
intersection = self.intersect(other)
# Check if the intersection is non-empty and has dimension (d-1)
if intersection.is_empty or intersection.affine_dimension != self.dim - 1:
return False
# Check if the intersection is a facet of both polyhedra
for i in range(self.A.shape[0]):
facet = self.get_facet(i)
if intersection in facet:
for j in range(other.A.shape[0]):
other_facet = other.get_facet(j)
if intersection in other_facet:
return True
return False
[docs]
def hausdorff_distance(self, other: "Polyhedron") -> float:
r"""Compute the Hausdorff distance between two polyhedra.
For two polyhedra :math:`P` and :math:`Q`, the Hausdorff distance is defined as:
.. math::
d_H(P, Q) = \max\left\{ \max\limits_{x\in P} \min\limits_{y\in Q} \|x-y\|, \max\limits_{y\in Q} \min\limits_{x\in P} \|x-y\| \right\}.
.. note::
This method is computationally expensive, as it requires to compute the distances of all vertices to another set for both polyhedra.
"""
if not isinstance(other, Polyhedron):
raise TypeError("The other object must be a Polyhedron.")
if self.dim != other.dim:
raise mpterr.DimMismatchError(self.dim, other.dim)
d1 = max(self.distance(vertex_other) for vertex_other in other.V)
d2 = max(other.distance(vertex_self) for vertex_self in self.V)
return float(max(d1, d2))
# Functions to implement
# Requires parametric solver
# - projection
# - affineMap - A * P + b
# - invAffineMap
# - str and repr
# - chebycenter
# - getFacet
# - outerApprox
# - shoot (max alpha s.t. x + alpha * d in Set)
# - slice
# - triangulate
# - volume
# - Minkowski difference
# - Minkowski sum
# - Set difference (requires Polyhedral arrays)
# - isFullDim
# Functions for later
# - meshGrid
# Functions to skip
# - facetInteriorPoints
# - interiorPoint
# - minAffineRep
# TODO: provide some standard polyhedra constructors, e.g., box, simplex, hexagon, etc.
[docs]
class Hexagon(Polyhedron):
"""
A class to represent a regular hexagon in 2D.
"""
def __init__(self, center: Vector = np.zeros(2), radius: float = 1.0, rotation: float = 0.0):
"""
Initialize a regular hexagon.
Args:
center: Center of the hexagon.
radius: Distance from the center to each vertex.
rotation: Counter-clockwise rotation angle (around the center) in radians.
"""
# rot_matrix = np.array([[np.cos(rotation), -np.sin(rotation)],
# [np.sin(rotation), np.cos(rotation)]])
angles = np.linspace(0, 2 * np.pi, 6, endpoint=False)
angles += rotation
vertices = np.array([[center[0] + radius * np.cos(angle),
center[1] + radius * np.sin(angle)] for angle in angles])
super().__init__(V=VData(V=vertices))
[docs]
class Simplex(Polyhedron):
r"""
Represent a standard :math:`k`-simplex.
"""
def __init__(self, k: int):
if k < 0:
raise ValueError("Order of a simplex must be at least 0.")
dim = k + 1
vertices = np.eye(dim) # standard basis vectors in R^(k+1)
A = -np.eye(dim)
b = np.zeros(dim)
Ae = np.ones((1, dim))
be = np.ones(1)
super().__init__(H=HData(A, b, Ae, be, is_minimal=True, preprocess=False),
V=VData(V=vertices, is_minimal=True, preprocess=False))
self._affine_dim = k
self._chebyshev_center = 1 / dim * np.ones(self.dim)
self._chebyshev_radius = 1 / dim
@property
def is_bounded(self):
return True
class RandomPolyhedronGenerator:
"""
A class to generate random polyhedra in H-representation or V-representation.
"""
@staticmethod
def random_hrep(dim: int, num_ineq: int, num_eq: int = 0, bounded: bool = True, nonempty: bool = True, seed: Optional[int] = None) -> Polyhedron:
"""
Generate a random polyhedron in H-representation.
Args:
dim: Dimension of the polyhedron.
num_ineq: Number of inequality constraints.
num_eq: Number of equality constraints.
bounded: If True, ensures the polyhedron is bounded.
nonempty: If True, ensures the polyhedron is non-empty.
seed: Random seed for reproducibility.
Returns:
A Polyhedron object in H-representation.
"""
if seed is not None:
np.random.seed(seed)
A = np.random.randn(num_ineq, dim)
b = np.random.rand(num_ineq) * 10 # Random positive bounds
Ae = np.random.randn(num_eq, dim) if num_eq > 0 else np.empty((0, dim))
be = np.random.rand(num_eq) * 10 if num_eq > 0 else np.empty((0,))
if bounded:
# Add box constraints to ensure boundedness
box_A = np.vstack([np.eye(dim), -np.eye(dim)])
box_b = np.ones(2 * dim) * 10
A = np.vstack([A, box_A])
b = np.hstack([b, box_b])
if nonempty:
# Ensure the polyhedron is non-empty by adding a feasible point
feasible_point = np.random.rand(dim) * 5 # Random feasible point
A = np.vstack([A, -np.eye(dim)])
b = np.hstack([b, -(-feasible_point)])
return Polyhedron(H=HData(A=A, b=b, Ae=Ae, be=be))
@staticmethod
def random_vrep(dim: int, num_vertices: int, num_rays: int = 0, seed: Optional[int] = None) -> Polyhedron:
"""
Generate a random polyhedron in V-representation.
Args:
dim: Dimension of the polyhedron.
num_vertices: Number of vertices.
num_rays: Number of rays.
bounded: If True, ensures the polyhedron is bounded (no rays).
seed: Random seed for reproducibility.
Returns:
A Polyhedron object in V-representation.
"""
if seed is not None:
np.random.seed(seed)
V = np.random.rand(num_vertices, dim) * 10 # Random vertices
R = np.random.randn(num_rays, dim)
return Polyhedron(V=VData(V=V, R=R))