import numpy as np
from dataclasses import dataclass
from numpy.typing import NDArray
from typing import Optional, Union
from scipy.spatial import distance_matrix
from .tolerances import tolerances
Scalar = Union[int, float, np.integer, np.floating]
Vector = NDArray[np.float64]
Matrix = NDArray[np.float64]
Bool = np.bool_ | bool
IntVector = NDArray[np.int64]
IntMatrix = NDArray[np.int64]
[docs]
@dataclass
class HData:
r"""The halfspace representation (H-representation) of a polyhedron, defined as:
.. math::
\{ x\in\mathbb{R}^n \mid A x \leq b, A_e x = b_e \}.
.. Note::
If ``preprocess`` is set as `True`, the ``__init__`` method performs pre-processing on the provided arguments.
This might cause the `A`, `b`, `Ae`, and `be` attributes to be different from what the user provided in arguments.
"""
A: Matrix
b: Vector
Ae: Matrix
be: Vector
# True if the representation is minimal
is_minimal: bool = False
def __init__(self, A: Optional[Matrix] = None, b: Optional[Vector] = None,
Ae: Optional[Matrix] = None, be: Optional[Vector] = None,
is_minimal: bool = False, preprocess: bool = True):
A = None if A is None else np.asarray(A)
b = None if b is None else np.asarray(b)
Ae = None if Ae is None else np.asarray(Ae)
be = None if be is None else np.asarray(be)
# Make sure the data is numeric ndarray (detect the case when data can be converted to ndarray but the dtype is object)
if not (A is None or np.issubdtype(A.dtype, np.number)):
raise TypeError("A must be convertable to numeric numpy ndarray.")
if not (b is None or np.issubdtype(b.dtype, np.number)):
raise TypeError("b must be convertable to numeric numpy ndarray.")
if not (Ae is None or np.issubdtype(Ae.dtype, np.number)):
raise TypeError("Ae must be convertable to numeric numpy ndarray.")
if not (be is None or np.issubdtype(be.dtype, np.number)):
raise TypeError("be must be convertable to numeric numpy ndarray.")
# Must either provide (A, b) or (Ae, be), or both
if not ( (A is not None and b is not None) or (Ae is not None and be is not None) ):
raise ValueError("Either (A, b) or (Ae, be) must be provided!")
# Either A and b are both provided, or they are both not provided
if (A is None) != (b is None):
raise ValueError("A and b must be provided together or neither should be provided.")
self.A = np.zeros((0, Ae.shape[1])) if A is None else A
self.b = np.zeros(0) if b is None else b
# Detect NaN in arguments
if np.isnan(self.A).any():
raise ValueError("A can not contain nan values.")
if np.isnan(self.b).any():
raise ValueError("b can not contain nan values.")
if self.A.ndim != 2:
raise ValueError("A must be a matrix")
if self.b.ndim != 1:
raise ValueError("b must be a vector")
if self.A.shape[0] != self.b.shape[0]:
raise ValueError("Incompatible dimensions: number of rows of A must match the number of elements in b")
# Either Ae and be are both provided, or they are both not provided
if (Ae is None) != (be is None):
raise ValueError("Ae and be must be provided together or neither should be provided.")
self.Ae = np.zeros((0, self.A.shape[1])) if Ae is None else Ae
self.be = np.zeros(0) if be is None else be
if np.isnan(self.Ae).any():
raise ValueError("Ae can not contain nan values.")
if np.isnan(self.be).any():
raise ValueError("be can not contain nan values.")
if self.Ae.ndim != 2:
raise ValueError("Ae must be a matrix")
if self.be.ndim != 1:
raise ValueError("be must be a vector")
if self.Ae.shape[0] != self.be.shape[0]:
raise ValueError("Incompatible dimensions: number of rows of Ae must match the number of elements in be")
if self.A.shape[1] != self.Ae.shape[1]:
raise ValueError("Incompatible dimensions: number of columns of Ae must match the number of columns in A")
self.is_minimal = is_minimal
if preprocess:
self.preprocess()
[docs]
def preprocess(self):
r"""
Preprocess the H-representation by removing redundant constraints and normalizing the constraints.
- Deal with zero rows in :math:`A`:
If at least one row :math:`a_i` in :math:`A` is almost zero and its corresponding :math:`b_i`
is strictly negative, the provided H-presentation will be treated as an empty set,
represented by :math:`\{ x\in\mathbb{R}^n \mid 0_{1\times n}^\top x \leq -1\}`.
If at least one row :math:`a_i` in :math:`A` is almost zero and its corresponding :math:`b_i`
is non-negative, this row will be treated as redundant and removed, unless it is the
only row in :math:`A`. In this case, the provided H-presentation will be treated as the full space,
represented by :math:`\{ x\in\mathbb{R}^n \mid 0_{1\times n}^\top x \leq 1\}`.
- Deal with zero rows in :math:`A_e`:
If at least one row :math:`a_{e,i}` in :math:`A_e` is almost zero and its corresponding :math:`b_{e,i}`
is non-zero, the provided H-presentation will be treated as an empty set,
represented by :math:`\{ x\in\mathbb{R}^n \mid 0_{1\times n}^\top x \leq -1\}`.
If at least one row :math:`a_{e,i}` in :math:`A_e` is almost zero and its corresponding :math:`b_{e,i}`
is zero, this row will be treated as redundant and removed.
- Deal with rows that are multiples of other rows in :math:`A` and :math:`A_e`.
For :math:`A`, if there exist :math:`i<j`, s.t. :math:`a_i = t a_j, ~b_i = t b_j` where :math:`t \neq 0`,
then :math:`a_j` is treated as redundant and will be removed. Similar for :math:`A_e`.
.. Note::
This method modifies the H-representation in place.
"""
# ------------------ Deal with zeros rows in Ae --------------------
zero_rows_Ae = np.linalg.norm(self.Ae, axis=1) <= tolerances['zero_normal']
nonzero_rows_be = np.abs(self.be) >= tolerances['zero']
zero_rows_be = ~nonzero_rows_be
# If at least one equality is like [0, 0] * [x; y] = alpha, where alpha is non-zero.
# Then it implies an empty set.
if np.any(zero_rows_Ae & nonzero_rows_be):
self.A = np.zeros((1, self.dim))
self.b = np.array([-1])
self.Ae = np.zeros((0, self.dim))
self.be = np.zeros((0,))
self.is_minimal = True
return
# If at least one equality is like [0, 0] * [x; y] = 0, then this equality always holds and can be removed.
if np.any(zero_rows_Ae & zero_rows_be):
mask_Ae_be = ~(zero_rows_Ae & zero_rows_be)
self.Ae = self.Ae[mask_Ae_be]
self.be = self.be[mask_Ae_be]
# --------------------------------------
# ------------------ Deal with zeros rows in A --------------------
zero_rows_A = np.linalg.norm(self.A, axis=1) <= tolerances['zero_normal']
neg_rows_b = self.b <= -tolerances['zero']
nonneg_rows_b = ~neg_rows_b
# If at least one inequality is like [0, 0] * [x; y] <= alpha, where alpha is strictly negative,
# then it implies an empty set. We only preserve A = [0, 0] and b = -1
if np.any(zero_rows_A & neg_rows_b):
self.A = np.zeros((1, self.dim))
self.b = np.array([-1])
self.Ae = np.zeros((0, self.dim))
self.be = np.zeros((0,))
self.is_minimal = True
return
# Now we can safely assume if there are zero rows in A, the corresponding b of these rows are non-negative.
# These constraints are like [0, 0] * [x; y] <= beta, where beta >=0.
# These inequalities always hold, so we remove these inequalities.
if np.any(zero_rows_A & nonneg_rows_b):
mask_A_b = ~(zero_rows_A & nonneg_rows_b)
self.A = self.A[mask_A_b]
self.b = self.b[mask_A_b]
# --------------------------------------
# remove the rows that are (nearly) multiple of others, e.g., x <= 1 and 2x <= 2
# an alternative is to use index_to_keep = np.unique(H_normalized, axis=0, return_index=True), but this does not allow adjust the tolerance.
H = np.column_stack((self.A, self.b))
H_normalized = H / np.linalg.norm(H, axis=1).reshape((-1,1))
H_dists = distance_matrix(H_normalized, H_normalized)
H_index_to_keep = np.ones(self.A.shape[0], dtype=bool)
for i in range(self.A.shape[0]):
for j in range(i+1, self.A.shape[0]):
if H_dists[i, j] < tolerances['zero']:
H_index_to_keep[j] = False
self.A, self.b = self.A[H_index_to_keep], self.b[H_index_to_keep]
He = np.column_stack((self.Ae, self.be))
He_normalized = He / np.linalg.norm(He, axis=1).reshape((-1,1))
He_dists = distance_matrix(He_normalized, He_normalized)
He_index_to_keep = np.ones(self.Ae.shape[0], dtype=bool)
for i in range(self.Ae.shape[0]):
for j in range(i+1, self.Ae.shape[0]):
if He_dists[i, j] < tolerances['zero']:
He_index_to_keep[j] = False
self.Ae, self.be = self.Ae[He_index_to_keep], self.be[He_index_to_keep]
# If after the above check, we get empty A and Ae, this means the H-rep is the full space R^n
# In this case set the H-rep as [0, 0, ..., 0] x <= 1
if self.A.size == 0 and self.Ae.size == 0:
self.A = np.zeros((1, self.dim))
self.b = np.array([1])
self.is_minimal = True
@property
def dim(self) -> int:
return self.A.shape[1]
@property
def is_cone(self):
"""
Returns True if the Polyhedron is a cone
"""
return np.all(np.abs(self.b) < tolerances['zero']) and np.all(np.abs(self.be) < tolerances['zero'])
@property
def homogenize(self) -> "HData":
A = np.hstack([self.A, -self.b[:, np.newaxis]])
b = np.zeros(self.A.shape[0])
Ae = np.hstack([self.Ae, -self.be[:, np.newaxis]])
be = np.zeros(self.Ae.shape[0])
return HData(A=A,b=b,Ae=Ae,be=be)
@property
def dehomogenize(self) -> "HData":
A = self.A[:, :-1]
b = -self.A[:, -1]
Ae = self.Ae[:, :-1]
be = -self.Ae[:, -1]
# TODO: if the homogenized H-rep is minimal, is the de-homogenized one also minimal? Vice versa?
return HData(A=A, b=b, Ae=Ae, be=be)
def __hash__(self):
r"""
Compute a hash value for the HData object. This is useful for using HData objects as keys in dictionaries
or adding them to sets. The hash is computed based on the contents of A, b, Ae, and be.
.. Note::
The hash function assumes that the arrays A, b, Ae, and be are immutable. Modifying these arrays
after the hash is computed may lead to unexpected behavior.
"""
A_hash = hash(self.A.tobytes())
b_hash = hash(self.b.tobytes())
Ae_hash = hash(self.Ae.tobytes())
be_hash = hash(self.be.tobytes())
is_minimal_hash = hash(self.is_minimal)
return hash((A_hash, b_hash, Ae_hash, be_hash, is_minimal_hash))
[docs]
@dataclass
class VData:
r"""The vertex representation (V-representation) of a polyhedron, defined as:
.. math::
\{ x\in\mathbb{R}^n \mid x = V^\top \lambda + R^\top \mu, \lambda,\mu \geq 0, 1^\top \lambda = 1 \}.
.. Note::
If ``preprocess`` is set as `True`, the ``__init__`` method performs pre-processing on the provided arguments.
This might cause the `V` and `R` attributes to be different from what the user provided in arguments.
"""
V: Matrix
R: Matrix
# True if the representation is minimal
is_minimal: bool = False
def __init__(self, V: Optional[Matrix] = None, R: Optional[Matrix] = None,
is_minimal: bool = False, preprocess: bool = True):
V = None if V is None else np.asarray(V)
R = None if R is None else np.asarray(R)
# Make sure the data is numeric ndarray (detect the case when data can be converted to ndarray but the dtype is object)
if not (V is None or np.issubdtype(V.dtype, np.number)):
raise TypeError("V must be convertable to numeric numpy ndarray.")
if not (R is None or np.issubdtype(R.dtype, np.number)):
raise TypeError("R must be convertable to numeric numpy ndarray.")
if V is not None and R is not None:
self.V = V
self.R = R
elif V is None and R is not None:
self.V = np.zeros((0, R.shape[1]))
self.R = R
elif R is None and V is not None:
self.V = V
self.R = np.zeros((0, V.shape[1]))
else:
raise ValueError("Either V or R must be given.")
self.is_minimal = is_minimal
if self.V.ndim != 2:
raise ValueError("V must be a matrix")
if self.R.ndim != 2:
raise ValueError("R must be a matrix")
if self.V.shape[1] != self.R.shape[1]:
raise ValueError("Incompatible dimensions: number of columns of V must match the number of columns of R")
if preprocess:
self.preprocess()
[docs]
def preprocess(self):
r"""
Preprocess the V-representation by removing obviously redundant vertices and rays.
- Deal with repeated rows in :math:`V`:
If there are repeated rows in :math:`V`, i.e., :math:`v_i = v_j` where :math:`i<j`, then :math:`v_j`
will be treated as a redundant vertex and removed.
- Deal with zero rows in :math:`R`:
If there are zero rows in :math:`R`, these zeros rows will be removed.
- Deal with rows that are positive multiples of others in :math:`R`:
If there exist :math:`i<j`, s.t. :math:`r_i = t r_j` where :math:`t > 0`,
:math:`r_j` will be treated as a redundant ray and removed.
"""
# Remove the repeated elements in vertices
V_dists = distance_matrix(self.V, self.V) # compute distance matrix
index_triu = np.triu_indices(V_dists.shape[0], k=1) # only upper triangular part needed
redundant = index_triu[1][V_dists[index_triu] < tolerances['zero']]
redundant = np.unique(redundant)
self.V = np.delete(self.V, redundant, axis=0)
# Remove the zero rows in rays
mask_R = np.linalg.norm(self.R, axis=1) > tolerances['zero_normal']
self.R = self.R[mask_R]
# Remove the rays that're (positive) multiple of others
R_normalized = self.R / np.linalg.norm(self.R, axis=1, keepdims=True)
R_dists = distance_matrix(R_normalized, R_normalized) # compute distance matrix
index_triu = np.triu_indices(R_dists.shape[0], k=1) # only upper triangular part needed
redundant = index_triu[1][R_dists[index_triu] < tolerances['zero']]
redundant = np.unique(redundant)
self.R = np.delete(self.R, redundant, axis=0)
@property
def dim(self) -> int:
return self.V.shape[1]
@property
def is_cone(self):
"""
Returns True if the Polyhedron is a cone
"""
return self.V.shape[0] == 1 and np.all(np.abs(self.V) < tolerances['zero'])
@property
def homogenize(self) -> "VData":
R = np.block([[self.V, np.ones((self.V.shape[0], 1))],
[self.R, np.zeros((self.R.shape[0], 1))]])
V = np.zeros((1, self.dim + 1))
return VData(V=V,R=R)
@property
def dehomogenize(self) -> "VData":
r"""
Be careful when calling this function. It only makes sense if the V-representation is homogenous.
"""
if self.V.size != 0 and not np.allclose(self.V, 0):
raise ValueError("The V-Representation has vertices so it is non-homogeneous, therefore cannot be dehomogenized.")
vertices = np.zeros((0, self.dim-1))
rays = np.zeros((0, self.dim-1))
for i in range(self.R.shape[0]):
ray = self.R[i, :]
if np.allclose(ray[-1], 1, tolerances['zero']):
vertices = np.vstack((vertices, np.delete(ray, -1)))
# rays = np.zeros((0, self.dim))
elif np.allclose(ray[-1], 0, tolerances['zero']):
rays = np.vstack((rays, np.delete(ray, -1)))
else:
raise ValueError("At least one element in the rays is neither 0 not 1. "
"The V-representation is non-homogeneous therefore cannot be dehomogenized.")
return VData(V=vertices, R=rays)
def __hash__(self):
r"""
Compute a hash value for the VData object. This is useful for using VData objects as keys in dictionaries
or adding them to sets. The hash is computed based on the contents of V, R, and is_minimal.
.. Note::
The hash function assumes that the arrays V and R are immutable. Modifying these arrays
after the hash is computed may lead to unexpected behavior.
"""
V_hash = hash(self.V.tobytes())
R_hash = hash(self.R.tobytes())
is_minimal_hash = hash(self.is_minimal)
return hash((V_hash, R_hash, is_minimal_hash))