Pular para conteúdo

tucoopy.geometry.polyhedron

Polyhedral sets (H-representation) for geometry modules.

This module defines PolyhedralSet, a small wrapper around linear constraints in the form:

  • inequalities: \(A_{ub} x \leq b_{ub}\)
  • equalities: \(A_{eq} x = b_{eq}\)
  • variable bounds: \(l_i \leq x_i \leq u_i\)

It provides core operations used across tucoopy.geometry:

  • feasibility checks and extracting a sample point (LP),
  • membership tests,
  • Chebyshev center (LP),
  • (small-\(n\)) vertex enumeration and projection.
Notes

This class is intentionally lightweight and relies on the configured LP backend. For large-dimensional polytopes, vertex enumeration is not feasible; helpers in this module are explicitly guarded by max_dim limits.

Examples:

A 2D triangle (simplex): x>=0, y>=0, x+y<=1.

>>> from tucoopy.geometry import PolyhedralSet
>>> P = PolyhedralSet.from_hrep(
...     A_ub=[[-1, 0], [0, -1], [1, 1]],
...     b_ub=[0, 0, 1],
... )
>>> P.extreme_points()
[[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]]

PolyhedralSet dataclass

Polyhedral set in H-representation with optional bounds.

Definition

A polyhedral set (polyhedron) in \(\mathbb{R}\) is described by linear equalities, linear inequalities, and coordinate bounds:

\[ A_{ub} x \le b_{ub}, \\ A_{eq} x = b_{eq}, \\ l_i \le x_i \le u_i. \]

This class is a lightweight geometry utility intended to support set-valued solution objects (core, epsilon-core, least-core, etc.) and visualization tasks such as: - feasibility checks, - sampling a feasible point, - Chebyshev center, - enumerating extreme points in small dimension (brute-force), - projecting extreme points to selected coordinates.

Notes
  • LP-based methods require an LP backend (SciPy/HiGHS or PuLP via your linprog_solve wrapper).
  • The extreme_points enumerator is exponential and intended only for small dimension (visualization).

Attributes:

Name Type Description
A_ub, b_ub

Inequality constraints A_ub x <= b_ub.

A_eq, b_eq

Equality constraints A_eq x = b_eq.

bounds list[tuple[float | None, float | None]]

Per-coordinate bounds (lb, ub); each side may be None.

Examples:

A 2D triangle: x>=0, y>=0, x+y<=1.

>>> P = PolyhedralSet.from_hrep(
...     A_ub=[[ -1,  0], [ 0, -1], [1, 1]],
...     b_ub=[  0,  0, 1],
... )
>>> P.contains([0.2, 0.3])
True
>>> P.contains([0.8, 0.3])
False

n_vars property

n_vars

Dimension (number of variables) inferred from bounds or constraints.

Returns:

Type Description
int

The inferred dimension n.

Raises:

Type Description
InvalidParameterError

If the dimension cannot be inferred (no bounds and empty constraint matrices).

Examples:

>>> from tucoopy.geometry import PolyhedralSet
>>> P = PolyhedralSet.from_hrep(bounds=[(0.0, 1.0), (None, None)])
>>> P.n_vars
2

affine_dimension

affine_dimension(*, tol=DEFAULT_GEOMETRY_TOL)

Estimate the affine dimension of the equality/bound hull.

This estimates the dimension of the affine subspace defined by:

  • equalities: \(A_{eq} x = b_{eq}\), and
  • fixed bounds: \(x_i = l_i = u_i\).

The returned value is:

\[\dim \approx n - \mathrm{rank}(M),\]

where \(M\) stacks the equality rows and the unit vectors for fixed coordinates.

Parameters:

Name Type Description Default
tol float

Numerical tolerance used in the rank estimate.

DEFAULT_GEOMETRY_TOL

Returns:

Type Description
int

Estimated affine dimension in \(\{0,\dots,n\}\).

Notes
  • This method does not check feasibility. For an empty set, the result is a diagnostic estimate for the constraint hull rather than the set itself.
  • If there are no equalities and no fixed bounds, this returns \(n\).

Examples:

A line segment in 2D has affine dimension 1:

>>> from tucoopy.geometry import PolyhedralSet
>>> P = PolyhedralSet.from_hrep(A_eq=[[1.0, 1.0]], b_eq=[1.0], bounds=[(0.0, 1.0), (0.0, 1.0)])
>>> P.affine_dimension()
1

chebyshev_center

chebyshev_center()

Chebyshev center via LP (maximum inscribed ball).

The Chebyshev center is the point maximizing the radius of an Euclidean ball contained in the polyhedron. For inequalities \(a_k^T x \leq b_k\), the ball constraint becomes:

\[ a_k^T x + \|a_k\|_2 \, r \le b_k, \qquad r \geq 0. \]

Equalities are kept as \(a^T x = b\) (they do not involve \(r\)).

Returns:

Type Description
(center, radius) | None

The center x and radius r of the largest inscribed ball, or None if infeasible.

Notes
  • If there are no inequalities (only equalities/bounds), a meaningful inscribed-ball radius is not defined in the same way. In that case we return any feasible point and \(r = +\infty\) as a sentinel.

Examples:

>>> P = PolyhedralSet.from_hrep(A_ub=[[1, 0], [0, 1], [-1, 0], [0, -1]], b_ub=[1, 1, 0, 0])
>>> center, r = P.chebyshev_center()
>>> r >= 0
True

chebyshev_center_with_diag

chebyshev_center_with_diag()

Chebyshev center via LP with diagnostics.

Returns:

Type Description
(center_radius, result)

center_radius is (x, r) or None if infeasible. result is the LP backend result object.

Examples:

>>> from tucoopy.geometry import PolyhedralSet
>>> P = PolyhedralSet.from_hrep(bounds=[(0.0, 1.0)])
>>> (cc, res) = P.chebyshev_center_with_diag()
>>> cc is None or len(cc[0]) == 1
True

contains

contains(x, *, tol=DEFAULT_GEOMETRY_TOL)

Check membership with tolerance.

Parameters:

Name Type Description Default
x Sequence[float]

Candidate point.

required
tol float

Tolerance for bound checks, inequalities, and equalities.

DEFAULT_GEOMETRY_TOL

Returns:

Type Description
bool

True if x satisfies all constraints within tolerance.

Notes
  • Inequalities are accepted if \(A_{ub} x \leq b_{ub} + \text{tol}\).
  • Equalities are accepted if \(|A_{eq} x - b_{eq}| \leq \text{tol}\) row-wise.

Examples:

>>> from tucoopy.geometry import PolyhedralSet
>>> P = PolyhedralSet.from_hrep(A_ub=[[1, 1]], b_ub=[1], bounds=[(0.0, None), (0.0, None)])
>>> P.contains([0.25, 0.25])
True
>>> P.contains([0.75, 0.75])
False

extreme_points

extreme_points(
    *,
    tol=DEFAULT_GEOMETRY_TOL,
    max_dim=DEFAULT_GEOMETRY_MAX_DIM,
)

Enumerate extreme points (vertices) for small-dimensional polytopes.

Method

A point is a vertex if it is the unique intersection of \(n\) linearly independent active constraints (equalities plus tight inequalities/bounds).

This brute-force enumerator proceeds as follows:

  1. Convert variable bounds to additional inequalities.
  2. Choose \(n - m_{eq}\) inequalities to add to the \(m_eq\) equalities, yielding a square linear system \(A x = b\).
  3. Solve this system.
  4. Keep the solution if it satisfies the original polyhedron constraints.
  5. De-duplicate solutions by quantizing coordinates.

Parameters:

Name Type Description Default
tol float

Numerical tolerance used both in the linear solve and membership tests.

DEFAULT_GEOMETRY_TOL
max_dim int

Safety limit: only dimensions n <= max_dim are supported.

DEFAULT_GEOMETRY_MAX_DIM

Returns:

Type Description
list[list[float]]

List of vertices, sorted lexicographically. Returns an empty list if no vertices are found.

Notes
  • This routine is exponential in the number of inequalities and is intended for visualization (e.g. \(n=2,3,4\)) rather than serious polyhedral computation.
  • Degenerate polytopes (many redundant/tied constraints) may require a larger tolerance for stable de-duplication.

Examples:

Triangle in 2D:

>>> P = PolyhedralSet.from_hrep(A_ub=[[-1,0],[0,-1],[1,1]], b_ub=[0,0,1])
>>> P.extreme_points()
[[0.0, 0.0], [0.0, 1.0], [1.0, 0.0]]

from_hrep classmethod

from_hrep(
    *,
    A_ub=None,
    b_ub=None,
    A_eq=None,
    b_eq=None,
    bounds=None,
)

Construct a polyhedron from an H-representation.

Parameters:

Name Type Description Default
A_ub Sequence[Sequence[float]] | None

Inequality constraints A_ub x <= b_ub. If omitted, no inequalities.

None
b_ub Sequence[Sequence[float]] | None

Inequality constraints A_ub x <= b_ub. If omitted, no inequalities.

None
A_eq Sequence[Sequence[float]] | None

Equality constraints A_eq x = b_eq. If omitted, no equalities.

None
b_eq Sequence[Sequence[float]] | None

Equality constraints A_eq x = b_eq. If omitted, no equalities.

None
bounds Sequence[tuple[float | None, float | None]] | None

Optional bounds (lb, ub) per coordinate.

None

Returns:

Type Description
PolyhedralSet

A normalized polyhedral set.

Raises:

Type Description
InvalidParameterError

If dimensions mismatch or the dimension cannot be inferred.

Examples:

>>> from tucoopy.geometry import PolyhedralSet
>>> # Unit square in 2D: 0<=x<=1, 0<=y<=1.
>>> P = PolyhedralSet.from_hrep(bounds=[(0.0, 1.0), (0.0, 1.0)])
>>> P.contains([0.5, 0.5])
True

is_bounded

is_bounded(*, tol=DEFAULT_GEOMETRY_TOL)

Check boundedness by optimizing each coordinate (LP-based).

A polyhedron is bounded iff every linear functional attains a finite maximum. As a practical diagnostic, this method checks whether each coordinate \(x_i\) has a finite minimum and maximum by solving \(2n\) LPs.

Parameters:

Name Type Description Default
tol float

Tolerance used to detect fixed bounds and to interpret solver output.

DEFAULT_GEOMETRY_TOL

Returns:

Type Description
bool

True if bounded (or empty), False if an unbounded direction is found.

Notes
  • Requires an LP backend at runtime.
  • For large n, this can be expensive (\(2n\) solves).
  • This method is intended mainly for guarding algorithms like hit-and-run.

Examples:

>>> from tucoopy.geometry import PolyhedralSet
>>> P = PolyhedralSet.from_hrep(bounds=[(0.0, 1.0), (0.0, 1.0)])
>>> P.is_bounded()
True

is_empty

is_empty()

Feasibility check via LP.

This solves a trivial feasibility LP (zero objective) and returns whether the LP solver finds any feasible point.

Returns:

Type Description
bool

True if the set is infeasible (empty); False otherwise.

Notes

Requires an LP backend at runtime (e.g. SciPy/HiGHS through linprog_solve).

Examples:

>>> from tucoopy.geometry import PolyhedralSet
>>> # Infeasible: x <= 0 and x >= 1.
>>> P = PolyhedralSet.from_hrep(A_ub=[[1], [-1]], b_ub=[0, -1])
>>> P.is_empty()
True

is_empty_with_diag

is_empty_with_diag()

Feasibility check via LP with diagnostics.

Returns:

Type Description
(is_empty, result)

Where result is the LP backend result object.

Examples:

>>> from tucoopy.geometry import PolyhedralSet
>>> P = PolyhedralSet.from_hrep(A_ub=[[1.0], [-1.0]], b_ub=[0.0, -1.0])
>>> empty, res = P.is_empty_with_diag()
>>> empty
True

project

project(
    dims,
    *,
    tol=DEFAULT_GEOMETRY_TOL,
    max_dim=DEFAULT_GEOMETRY_MAX_DIM,
    approx_n_points=None,
    approx_seed=None,
)

Project enumerated extreme points to selected coordinates (visualization helper).

Parameters:

Name Type Description Default
dims Sequence[int]

Indices of coordinates to keep.

required
tol float

Passed to extreme_points.

DEFAULT_GEOMETRY_TOL
max_dim int

Passed to extreme_points.

DEFAULT_GEOMETRY_MAX_DIM
approx_n_points int | None

If provided and n_vars > max_dim, use hit-and-run sampling to produce approximately projected points. This requires boundedness and an LP backend to find a starting point.

None
approx_seed int | None

Optional RNG seed used by hit-and-run when approx_n_points is set.

None

Returns:

Type Description
list[list[float]]

Projected vertex list.

Notes

This is not a true polyhedral projection (which would require eliminating variables). It simply enumerates vertices in the full space (small \(n\)) and returns the selected coordinates.

Examples:

>>> from tucoopy.geometry import PolyhedralSet
>>> P = PolyhedralSet.from_hrep(bounds=[(0.0, 1.0), (0.0, 1.0), (0.0, 1.0)])
>>> # Project cube vertices to the first two coordinates.
>>> verts2 = P.project((0, 1), max_dim=3)
>>> [0.0, 0.0] in verts2 and [1.0, 1.0] in verts2
True

residual_eq

residual_eq(x)

Compute equality residuals \(a_k^T x - b_k\) for each equality row.

Parameters:

Name Type Description Default
x Sequence[float]

Point in R^n.

required

Returns:

Type Description
list[float]

Residuals for each equality constraint.

Examples:

>>> from tucoopy.geometry import PolyhedralSet
>>> P = PolyhedralSet.from_hrep(A_eq=[[1, 1]], b_eq=[1], bounds=[(0.0, 1.0), (0.0, 1.0)])
>>> P.residual_eq([0.25, 0.75])
[0.0]

sample_point

sample_point()

Find any feasible point via LP.

Returns:

Type Description
list[float] | None

A feasible point if one exists, otherwise None.

Notes

This is a convenience wrapper around a feasibility LP with zero objective.

Examples:

>>> from tucoopy.geometry import PolyhedralSet
>>> P = PolyhedralSet.from_hrep(bounds=[(0.0, 1.0), (0.0, 1.0)])
>>> x = P.sample_point()
>>> x is None or P.contains(x)
True

sample_point_with_diag

sample_point_with_diag()

Find any feasible point via LP with diagnostics.

Returns:

Type Description
(x, result)

x is a feasible point if one exists, otherwise None. result is the LP backend result object.

Examples:

>>> from tucoopy.geometry import PolyhedralSet
>>> P = PolyhedralSet.from_hrep(bounds=[(0.0, 1.0)])
>>> x, res = P.sample_point_with_diag()
>>> x is None or P.contains(x)
True

sample_points_hit_and_run

sample_points_hit_and_run(
    n_points,
    *,
    start_point=None,
    burn_in=DEFAULT_HIT_AND_RUN_BURN_IN,
    thinning=DEFAULT_HIT_AND_RUN_THINNING,
    seed=None,
    tol=DEFAULT_HIT_AND_RUN_TOL,
)

Sample points using hit-and-run within the polyhedron (approximate uniform sampling).

This sampler: 1) starts from a feasible point (either start_point or sample_point), 2) repeatedly samples a random direction (restricted to the equality subspace), 3) computes the feasible segment along that direction using all inequalities (including bounds encoded as inequalities), 4) samples uniformly on that segment.

Parameters:

Name Type Description Default
n_points int

Number of returned samples.

required
start_point Sequence[float] | None

Optional feasible start point. If omitted, a start point is found via sample_point (which requires an LP backend, e.g. tucoopy[lp]).

None
burn_in int

Number of initial steps discarded.

DEFAULT_HIT_AND_RUN_BURN_IN
thinning int

Keep one sample every thinning steps after burn-in.

DEFAULT_HIT_AND_RUN_THINNING
seed int | None

Random seed.

None
tol float

Numerical tolerance used in feasibility checks.

DEFAULT_HIT_AND_RUN_TOL

Returns:

Type Description
list[list[float]]

Sampled points (each a length-n vector).

Notes
  • Requires a feasible start point. If start_point is not provided, this method calls sample_point, which requires an LP backend.
  • If the polyhedron is unbounded in a sampled direction (infinite segment), this method raises ValueError because "uniform" sampling is not defined.

Examples:

>>> from tucoopy.geometry import PolyhedralSet
>>> P = PolyhedralSet.from_hrep(bounds=[(0.0, 1.0), (0.0, 1.0)])
>>> pts = P.sample_points_hit_and_run(3, seed=0)
>>> len(pts)
3

slack_eq

slack_eq(x)

Compute absolute equality slacks \(|a_k^T x - b_k|\) for each equality row.

Parameters:

Name Type Description Default
x Sequence[float]

Point in R^n.

required

Returns:

Type Description
list[float]

Absolute residuals for each equality constraint.

Examples:

>>> from tucoopy.geometry import PolyhedralSet
>>> P = PolyhedralSet.from_hrep(A_eq=[[1.0, 1.0]], b_eq=[1.0])
>>> P.slack_eq([0.25, 0.75])
[0.0]

slack_ub

slack_ub(x)

Compute inequality slacks for the combined \(A_{ub} x \leq b_{ub}\) system.

This includes bound-inequalities (lb/ub converted to halfspaces) so it is suitable for "tight set" diagnostics.

Parameters:

Name Type Description Default
x Sequence[float]

Point in R^n.

required

Returns:

Type Description
list[float]

Slacks b_k - a_k^T x for each inequality row k.

Examples:

>>> from tucoopy.geometry import PolyhedralSet
>>> P = PolyhedralSet.from_hrep(A_ub=[[1]], b_ub=[1], bounds=[(0.0, None)])
>>> P.slack_ub([0.25])[0]
0.75