Common Set Representations

Ball2{N<:AbstractFloat} <: AbstractCentrallySymmetric{N}

Type that represents a ball in the 2-norm.

Fields

• center – center of the ball as a real vector
• radius – radius of the ball as a real scalar ($≥ 0$)

Notes

Mathematically, a ball in the 2-norm is defined as the set

where $c ∈ \mathbb{R}^n$ is its center and $r ∈ \mathbb{R}_+$ its radius. Here $‖ ⋅ ‖_2$ denotes the Euclidean norm (also known as 2-norm), defined as $‖ x ‖_2 = \left( \sum\limits_{i=1}^n |x_i|^2 \right)^{1/2}$ for any $x ∈ \mathbb{R}^n$.

Examples

Create a five-dimensional ball B in the 2-norm centered at the origin with radius 0.5:

julia> B = Ball2(zeros(5), 0.5)
Ball2{Float64}([0.0, 0.0, 0.0, 0.0, 0.0], 0.5)
julia> dim(B)
5

Evaluate B's support vector in the direction $[1,2,3,4,5]$:

julia> σ([1.,2.,3.,4.,5.], B)
5-element Array{Float64,1}:
 0.06741998624632421
 0.13483997249264842
 0.20225995873897262
 0.26967994498529685
 0.3370999312316211

σ(d::AbstractVector{N}, B::Ball2{N}) where {N<:AbstractFloat}

Return the support vector of a 2-norm ball in a given direction.

Input

• d – direction
• B – ball in the 2-norm

Output

The support vector in the given direction. If the direction has norm zero, the origin is returned.

Notes

Let $c$ and $r$ be the center and radius of a ball $B$ in the 2-norm, respectively. For nonzero direction $d$ we have

This function requires computing the 2-norm of the input direction, which is performed in the given precision of the numeric datatype of both the direction and the set. Exact inputs are not supported.

∈(x::AbstractVector{N}, B::Ball2{N})::Bool where {N<:AbstractFloat}

Check whether a given point is contained in a ball in the 2-norm.

Input

• x – point/vector
• B – ball in the 2-norm

Output

true iff $x ∈ B$.

Notes

This implementation is worst-case optimized, i.e., it is optimistic and first computes (see below) the whole sum before comparing to the radius. In applications where the point is typically far away from the ball, a fail-fast implementation with interleaved comparisons could be more efficient.

Algorithm

Let $B$ be an $n$-dimensional ball in the 2-norm with radius $r$ and let $c_i$ and $x_i$ be the ball's center and the vector $x$ in dimension $i$, respectively. Then $x ∈ B$ iff $\left( ∑_{i=1}^n |c_i - x_i|^2 \right)^{1/2} ≤ r$.

Examples

julia> B = Ball2([1., 1.], sqrt(0.5))
Ball2{Float64}([1.0, 1.0], 0.7071067811865476)
julia> [.5, 1.6] ∈ B
false
julia> [.5, 1.5] ∈ B
true

center(B::Ball2{N})::Vector{N} where {N<:AbstractFloat}

Return the center of a ball in the 2-norm.

Input

• B – ball in the 2-norm

Output

The center of the ball in the 2-norm.

rand(::Type{Ball2}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::Ball2{N}

Create a random ball in the 2-norm.

Input

• Ball2 – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A random ball in the 2-norm.

Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1. Additionally, the radius is nonnegative.

sample(B::Ball2{N}, nsamples::Int=1;
       [rng]::AbstractRNG=GLOBAL_RNG,
       [seed]::Union{Int, Nothing}=nothing) where {N<:AbstractFloat}

Return samples from a uniform distribution on the given ball in the 2-norm.

Input

• B – ball in the 2-norm
• nsamples – (optional, default: 1) number of samples
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A linear array of nsamples elements drawn from a uniform distribution in B.

Algorithm

Random sampling with uniform distribution in B is computed using Muller's method of normalized Gaussians. This function requires the package Distributions. See _sample_unit_nball_muller! for implementation details.

translate(B::Ball2{N}, v::AbstractVector{N}) where {N<:AbstractFloat}

Translate (i.e., shift) a ball in the 2-norm by a given vector.

Input

• B – ball in the 2-norm
• v – translation vector

Output

A translated ball in the 2-norm.

Algorithm

We add the vector to the center of the ball.

BallInf{N<:Real} <: AbstractHyperrectangle{N}

Type that represents a ball in the infinity norm.

Fields

• center – center of the ball as a real vector
• radius – radius of the ball as a real scalar ($≥ 0$)

Notes

Mathematically, a ball in the infinity norm is defined as the set

where $c ∈ \mathbb{R}^n$ is its center and $r ∈ \mathbb{R}_+$ its radius. Here $‖ ⋅ ‖_∞$ denotes the infinity norm, defined as $‖ x ‖_∞ = \max\limits_{i=1,…,n} \vert x_i \vert$ for any $x ∈ \mathbb{R}^n$.

Examples

Create the two-dimensional unit ball and compute its support function along the positive $x=y$ direction:

julia> B = BallInf(zeros(2), 1.0)
BallInf{Float64}([0.0, 0.0], 1.0)
julia> dim(B)
2
julia> ρ([1., 1.], B)
2.0

center(B::BallInf{N})::Vector{N} where {N<:Real}

Return the center of a ball in the infinity norm.

Input

• B – ball in the infinity norm

Output

The center of the ball in the infinity norm.

radius(B::BallInf, [p]::Real=Inf)::Real

Return the radius of a ball in the infinity norm.

Input

• B – ball in the infinity norm
• p – (optional, default: Inf) norm

Output

A real number representing the radius.

Notes

The radius is defined as the radius of the enclosing ball of the given $p$-norm of minimal volume with the same center.

radius_hyperrectangle(B::BallInf{N})::Vector{N} where {N<:Real}

Return the box radius of a ball in the infinity norm, which is the same in every dimension.

Input

• B – ball in the infinity norm

Output

The box radius of the ball in the infinity norm.

radius_hyperrectangle(B::BallInf{N}, i::Int)::N where {N<:Real}

Return the box radius of a ball in the infinity norm in a given dimension.

Input

• B – ball in the infinity norm
• i – dimension of interest

Output

The box radius of the ball in the infinity norm in the given dimension.

isflat(B::BallInf)::Bool

Determine whether a ball in the infinity norm is flat, i.e. whether its radius is zero.

Input

• B – ball in the infinity norm

Output

true iff the ball is flat.

Notes

For robustness with respect to floating-point inputs, this function relies on the result of isapproxzero when applied to the radius of the ball. Hence, this function depends on the absolute zero tolerance ABSZTOL.

rand(::Type{BallInf}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::BallInf{N}

Create a random ball in the infinity norm.

Input

• BallInf – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A random ball in the infinity norm.

Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1. Additionally, the radius is nonnegative.

σ(d::AbstractVector{N}, B::BallInf{N}) where {N<:Real}

Return the support vector of a ball in the infinity norm in a given direction.

Input

• d – direction
• B – ball in the infinity norm

Output

The support vector in the given direction. If the direction has norm zero, the vertex with biggest values is returned.

ρ(d::AbstractVector{N}, B::BallInf{N}) where {N<:Real}

Evaluate the support function of a ball in the infinity norm in a given direction.

Input

• d – direction
• B – ball in the infinity norm

Output

Evaluation of the support function in the given direction.

Algorithm

Let $B$ be a ball in the infinity norm with center $c$ and radius $r$ and let d be the direction of interest. For balls with dimensions less than 30 we use the implementation for AbstractHyperrectangle, taylored to a BallInf, which computes

where $\sgn(α) = 1$ if $α ≥ 0$ and $\sgn(α) = 1$ if $α < 0$.

For balls of higher dimension, we instead exploit that for a support vector $v = σ(d, B) = c + \sgn(d) * (r, …, r)ᵀ$ we have

where $⟨·, ·⟩$ denotes the dot product.

translate(B::BallInf{N}, v::AbstractVector{N}) where {N<:Real}

Translate (i.e., shift) a ball in the infinity norm by a given vector.

Input

• B – ball in the infinity norm
• v – translation vector

Output

A translated ball in the infinity norm.

Algorithm

We add the vector to the center of the ball.

Ball1{N<:Real} <: AbstractCentrallySymmetricPolytope{N}

Type that represents a ball in the 1-norm (also known as the Manhattan norm). The ball is also known as a cross-polytope.

It is defined as the set

where $c ∈ \mathbb{R}^n$ is its center and $r ∈ \mathbb{R}_+$ its radius.

Fields

• center – center of the ball as a real vector
• radius – radius of the ball as a scalar ($≥ 0$)

Examples

Unit ball in the 1-norm in the plane:

julia> B = Ball1(zeros(2), 1.)
Ball1{Float64}([0.0, 0.0], 1.0)
julia> dim(B)
2

We evaluate the support vector in the East direction:

julia> σ([0.,1], B)
2-element Array{Float64,1}:
 0.0
 1.0

σ(d::AbstractVector{N}, B::Ball1{N}) where {N<:Real}

Return the support vector of a ball in the 1-norm in a given direction.

Input

• d – direction
• B – ball in the 1-norm

Output

Support vector in the given direction.

∈(x::AbstractVector{N}, B::Ball1{N})::Bool where {N<:Real}

Check whether a given point is contained in a ball in the 1-norm.

Input

• x – point/vector
• B – ball in the 1-norm

Output

true iff $x ∈ B$.

Notes

This implementation is worst-case optimized, i.e., it is optimistic and first computes (see below) the whole sum before comparing to the radius. In applications where the point is typically far away from the ball, a fail-fast implementation with interleaved comparisons could be more efficient.

Algorithm

Let $B$ be an $n$-dimensional ball in the 1-norm with radius $r$ and let $c_i$ and $x_i$ be the ball's center and the vector $x$ in dimension $i$, respectively. Then $x ∈ B$ iff $∑_{i=1}^n |c_i - x_i| ≤ r$.

Examples

julia> B = Ball1([1., 1.], 1.);

julia> [.5, -.5] ∈ B
false
julia> [.5, 1.5] ∈ B
true

vertices_list(B::Ball1{N})::Vector{Vector{N}} where {N<:Real}

Return the list of vertices of a ball in the 1-norm.

Input

• B – ball in the 1-norm

Output

A list containing the vertices of the ball in the 1-norm.

center(B::Ball1{N})::Vector{N} where {N<:Real}

Return the center of a ball in the 1-norm.

Input

• B – ball in the 1-norm

Output

The center of the ball in the 1-norm.

rand(::Type{Ball1}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::Ball1{N}

Create a random ball in the 1-norm.

Input

• Ball1 – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A random ball in the 1-norm.

Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1. Additionally, the radius is nonnegative.

constraints_list(P::Ball1{N}) where {N<:Real}

Return the list of constraints defining a ball in the 1-norm.

Input

• B – ball in the 1-norm

Output

The list of constraints of the ball.

Algorithm

The constraints can be defined as $d_i^T (x-c) ≤ r$ for all $d_i$, where $d_i$ is a vector with elements $1$ or $-1$ in $n$ dimensions. To span all possible $d_i$, the function Iterators.product is used.

translate(B::Ball1{N}, v::AbstractVector{N}) where {N<:Real}

Translate (i.e., shift) a ball in the 1-norm by a given vector.

Input

• B – ball in the 1-norm
• v – translation vector

Output

A translated ball in the 1-norm.

Algorithm

We add the vector to the center of the ball.

Ballp{N<:AbstractFloat} <: AbstractCentrallySymmetric{N}

Type that represents a ball in the p-norm, for $1 ≤ p ≤ ∞$.

It is defined as the set

where $c ∈ \mathbb{R}^n$ is its center and $r ∈ \mathbb{R}_+$ its radius. Here $‖ ⋅ ‖_p$ for $1 ≤ p ≤ ∞$ denotes the vector $p$-norm, defined as $‖ x ‖_p = \left( \sum\limits_{i=1}^n |x_i|^p \right)^{1/p}$ for any $x ∈ \mathbb{R}^n$.

Fields

• p – norm as a real scalar
• center – center of the ball as a real vector
• radius – radius of the ball as a scalar ($≥ 0$)

Notes

The special cases $p=1$, $p=2$ and $p=∞$ fall back to the specialized types Ball1, Ball2 and BallInf, respectively.

Examples

A five-dimensional ball in the $p=3/2$ norm centered at the origin of radius 0.5:

julia> B = Ballp(3/2, zeros(5), 0.5)
Ballp{Float64}(1.5, [0.0, 0.0, 0.0, 0.0, 0.0], 0.5)
julia> dim(B)
5

We evaluate the support vector in direction $[1,2,…,5]$:

julia> σ([1., 2, 3, 4, 5], B)
5-element Array{Float64,1}:
 0.013516004434607558
 0.05406401773843023
 0.12164403991146802
 0.21625607095372093
 0.33790011086518895

σ(d::AbstractVector{N}, B::Ballp{N}) where {N<:AbstractFloat}

Return the support vector of a Ballp in a given direction.

Input

• d – direction
• B – ball in the p-norm

Output

The support vector in the given direction. If the direction has norm zero, the center of the ball is returned.

Algorithm

The support vector of the unit ball in the $p$-norm along direction $d$ is:

where $\tilde{v}_i = \frac{|d_i|^q}{d_i}$ if $d_i ≠ 0$ and $tilde{v}_i = 0$ otherwise, for all $i=1,…,n$, and $q$ is the conjugate number of $p$. By the affine transformation $x = r\tilde{x} + c$, one obtains that the support vector of $\mathcal{B}_p^n(c, r)$ is

where $v_i = c_i + r\frac{|d_i|^q}{d_i}$ if $d_i ≠ 0$ and $v_i = 0$ otherwise, for all $i = 1, …, n$.

∈(x::AbstractVector{N}, B::Ballp{N})::Bool where {N<:AbstractFloat}

Check whether a given point is contained in a ball in the p-norm.

Input

• x – point/vector
• B – ball in the p-norm

Output

true iff $x ∈ B$.

Notes

This implementation is worst-case optimized, i.e., it is optimistic and first computes (see below) the whole sum before comparing to the radius. In applications where the point is typically far away from the ball, a fail-fast implementation with interleaved comparisons could be more efficient.

Algorithm

Let $B$ be an $n$-dimensional ball in the p-norm with radius $r$ and let $c_i$ and $x_i$ be the ball's center and the vector $x$ in dimension $i$, respectively. Then $x ∈ B$ iff $\left( ∑_{i=1}^n |c_i - x_i|^p \right)^{1/p} ≤ r$.

Examples

julia> B = Ballp(1.5, [1., 1.], 1.)
Ballp{Float64}(1.5, [1.0, 1.0], 1.0)
julia> [.5, -.5] ∈ B
false
julia> [.5, 1.5] ∈ B
true

center(B::Ballp{N})::Vector{N} where {N<:AbstractFloat}

Return the center of a ball in the p-norm.

Input

• B – ball in the p-norm

Output

The center of the ball in the p-norm.

rand(::Type{Ballp}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::Ballp{N}

Create a random ball in the p-norm.

Input

• Ballp – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A random ball in the p-norm.

Algorithm

The center and radius are normally distributed with mean 0 and standard deviation 1. Additionally, the radius is nonnegative. The p-norm is a normally distributed number ≥ 1 with mean 1 and standard deviation 1.

translate(B::Ballp{N}, v::AbstractVector{N}) where {N<:AbstractFloat}

Translate (i.e., shift) a ball in the p-norm by a given vector.

Input

• B – ball in the p-norm
• v – translation vector

Output

A translated ball in the p- norm.

Algorithm

We add the vector to the center of the ball.

Ellipsoid{N<:AbstractFloat} <:  AbstractCentrallySymmetric{N}

Type that represents an ellipsoid.

It is defined as the set

where $c \in \mathbb{R}^n$ is its center and $Q \in \mathbb{R}^{n×n}$ its shape matrix, which should be a positive definite matrix. An ellipsoid can also be characterized as the image of a Euclidean ball by an invertible linear transformation. It is the higher-dimensional generalization of an ellipse.

Fields

• center – center of the ellipsoid
• shape_matrix – real positive definite matrix, i.e. it is equal to its transpose and $x^\mathrm{T}Qx > 0$ for all nonzero $x$

Examples

If the center is not specified, it is assumed that the center is the origin. For instance, a 3D ellipsoid with center at the origin and the shape matrix being the identity can be created with:

julia> using LinearAlgebra

julia> E = Ellipsoid(Matrix{Float64}(I, 3, 3))
Ellipsoid{Float64}([0.0, 0.0, 0.0], [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0])

julia> dim(E)
3

julia> an_element(E)
3-element Array{Float64,1}:
 0.0
 0.0
 0.0

This ellipsoid corresponds to the unit Euclidean ball. Let's evaluate its support vector in a given direction:

julia> σ(ones(3), E)
3-element Array{Float64,1}:
 0.5773502691896258
 0.5773502691896258
 0.5773502691896258

A two-dimensional ellipsoid with given center and shape matrix:

julia> E = Ellipsoid(ones(2), Diagonal([2.0, 0.5]))
Ellipsoid{Float64}([1.0, 1.0], [2.0 0.0; 0.0 0.5])

ρ(d::AbstractVector{N}, E::Ellipsoid{N}) where {N<:AbstractFloat}

Return the support function of an ellipsoid in a given direction.

Input

• d – direction
• E – ellipsoid

Output

The support function of the ellipsoid in the given direction.

Algorithm

The support value is $cᵀ d + ‖Bᵀ d‖₂$ where $c$ is the center and $Q = B Bᵀ$ is the shape matrix of E.

σ(d::AbstractVector{N}, E::Ellipsoid{N}) where {N<:AbstractFloat}

Return the support vector of an ellipsoid in a given direction.

Input

• d – direction
• E – ellipsoid

Output

Support vector in the given direction.

Algorithm

Let $E$ be an ellipsoid of center $c$ and shape matrix $Q = BB^\mathrm{T}$. Its support vector along direction $d$ can be deduced from that of the unit Euclidean ball $\mathcal{B}_2$ using the algebraic relations for the support vector,

∈(x::AbstractVector{N}, E::Ellipsoid{N})::Bool where {N<:AbstractFloat}

Check whether a given point is contained in an ellipsoid.

Input

• x – point/vector
• E – ellipsoid

Output

true iff x ∈ E.

Algorithm

The point $x$ belongs to the ellipsoid of center $c$ and shape matrix $Q$ if and only if

rand(::Type{Ellipsoid}; [N]::Type{<:AbstractFloat}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::Ellipsoid{N}

Create a random ellipsoid.

Input

• Ellipsoid – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A random ellipsoid.

Algorithm

The center is a normally distributed vector with entries of mean 0 and standard deviation 1.

The idea for the shape matrix comes from here. The matrix is symmetric positive definite, but also diagonally dominant.

where $n$ = dim (defaults to 2), and $S$ is a $n \times n$ random matrix whose coefficients are uniformly distributed in the interval $[-1, 1]$.

center(E::Ellipsoid{N})::Vector{N} where {N<:AbstractFloat}

Return the center of the ellipsoid.

Input

• E – ellipsoid

Output

The center of the ellipsoid.

translate(E::Ellipsoid{N}, v::AbstractVector{N}; share::Bool=false
         ) where {N<:AbstractFloat}

Translate (i.e., shift) an ellipsoid by a given vector.

Input

• E – ellipsoid
• v – translation vector
• share – (optional, default: false) flag for sharing unmodified parts of the original set representation

Output

A translated ellipsoid.

Notes

The shape matrix is shared with the original ellipsoid if share == true.

Algorithm

We add the vector to the center of the ellipsoid.

EmptySet{N<:Real} <: LazySet{N}

Type that represents the empty set, i.e., the set with no elements.

∅

An EmptySet instance of type Float64.

dim(∅::EmptySet)

Return the dimension of the empty set, which is -1 by convention.

Input

• ∅ – an empty set

Output

-1 by convention.

σ(d::AbstractVector{N}, ∅::EmptySet{N}) where {N<:Real}

Return the support vector of an empty set.

Input

• d – direction
• ∅ – empty set

Output

An error.

ρ(d::AbstractVector{N}, ∅::EmptySet{N}) where {N<:Real}

Evaluate the support function of an empty set in a given direction.

Input

• d – direction
• ∅ – empty set

Output

An error.

∈(x::AbstractVector{N}, ∅::EmptySet{N})::Bool where {N<:Real}

Check whether a given point is contained in an empty set.

Input

• x – point/vector
• ∅ – empty set

Output

The output is always false.

Examples

julia> [1.0, 0.0] ∈ ∅
false

an_element(∅::EmptySet)

Return some element of an empty set.

Input

• ∅ – empty set

Output

An error.

rand(::Type{EmptySet}; [N]::Type{<:Real}=Float64, [dim]::Int=0,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::EmptySet{N}

Create an empty set (note that there is nothing to randomize).

Input

• EmptySet – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 0) dimension (is ignored)
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

The (only) empty set of the given numeric type.

isbounded(∅::EmptySet)::Bool

Determine whether an empty set is bounded.

Input

• ∅ – empty set

Output

true.

isempty(∅::EmptySet)::Bool

Return if the empty set is empty or not.

Input

• ∅ – empty set

Output

true.

norm(S::EmptySet, [p]::Real=Inf)

Return the norm of an empty set. It is the norm of the enclosing ball (of the given $p$-norm) of minimal volume that is centered in the origin.

Input

• S – empty set
• p – (optional, default: Inf) norm

Output

An error.

radius(S::EmptySet, [p]::Real=Inf)

Return the radius of an empty set. It is the radius of the enclosing ball (of the given $p$-norm) of minimal volume with the same center.

Input

• S – empty set
• p – (optional, default: Inf) norm

Output

An error.

diameter(S::EmptySet, [p]::Real=Inf)

Return the diameter of an empty set. It is the maximum distance between any two elements of the set, or, equivalently, the diameter of the enclosing ball (of the given $p$-norm) of minimal volume with the same center.

Input

• S – empty set
• p – (optional, default: Inf) norm

Output

An error.

linear_map(M::AbstractMatrix{N}, ∅::EmptySet{N}) where {N}

Return the linear map of an empty set.

Input

• M – matrix
• ∅ – empty set

Output

The empty set.

translate(∅::EmptySet{N}, v::AbstractVector{N}) where {N<:Real}

Translate (i.e., shift) an empty set by a given vector.

Input

• ∅ – empty set
• v – translation vector

Output

The empty set.

plot_recipe(∅::EmptySet{N}, [ε]::N=zero(N)) where {N<:Real}

Convert an empty set to a sequence of points for plotting. In the special case of an empty set, we define the sequence as nothing.

Input

• ∅ – empty set
• ε – (optional, default: 0) ignored, used for dispatch

Output

nothing.

plot_emptyset(∅::EmptySet, [ε]::N=zero(N); ...)

Plot an empty set.

Input

• ∅ – empty set
• ε – (optional, default: 0) ignored, used for dispatch

HalfSpace{N<:Real, VN<:AbstractVector{N}} <: AbstractPolyhedron{N}

Type that represents a (closed) half-space of the form $a⋅x ≤ b$.

Fields

• a – normal direction (non-zero)
• b – constraint

Examples

The set $y ≥ 0$ in the plane:

julia> HalfSpace([0, -1.], 0.)
HalfSpace{Float64,Array{Float64,1}}([0.0, -1.0], 0.0)

LinearConstraint

Alias for HalfSpace

dim(hs::HalfSpace)::Int

Return the dimension of a half-space.

Input

• hs – half-space

Output

The ambient dimension of the half-space.

ρ(d::AbstractVector{N}, hs::HalfSpace{N})::N where {N<:Real}

Evaluate the support function of a half-space in a given direction.

Input

• d – direction
• hs – half-space

Output

The support function of the half-space. If the set is unbounded in the given direction, the result is Inf.

σ(d::AbstractVector{N}, hs::HalfSpace{N}) where {N<:Real}

Return the support vector of a half-space.

Input

• d – direction
• hs – half-space

Output

The support vector in the given direction, which is only defined in the following two cases:

1. The direction has norm zero.
2. The direction is the half-space's normal direction.

In both cases the result is any point on the boundary (the defining hyperplane). Otherwise this function throws an error.

∈(x::AbstractVector{N}, hs::HalfSpace{N})::Bool where {N<:Real}

Check whether a given point is contained in a half-space.

Input

• x – point/vector
• hs – half-space

Output

true iff $x ∈ hs$.

Algorithm

We just check if $x$ satisfies $a⋅x ≤ b$.

an_element(hs::HalfSpace{N})::Vector{N} where {N<:Real}

Return some element of a half-space.

Input

• hs – half-space

Output

An element on the defining hyperplane.

rand(::Type{HalfSpace}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::HalfSpace{N}

Create a random half-space.

Input

• HalfSpace – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A random half-space.

Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1. Additionally, the constraint a is nonzero.

isbounded(hs::HalfSpace)::Bool

Determine whether a half-space is bounded.

Input

• hs – half-space

Output

false.

isuniversal(hs::HalfSpace{N}, [witness]::Bool=false
           )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether a half-space is universal.

Input

• P – half-space
• witness – (optional, default: false) compute a witness if activated

Output

• If witness option is deactivated: false
• If witness option is activated: (false, v) where $v ∉ P$

Algorithm

Witness production falls back to isuniversal(::Hyperplane).

isempty(hs::HalfSpace)::Bool

Return if a half-space is empty or not.

Input

• hs – half-space

Output

false.

constraints_list(hs::HalfSpace{N}) where {N<:Real}

Return the list of constraints of a half-space.

Input

• hs – half-space

Output

A singleton list containing the half-space.

constraints_list(A::AbstractMatrix{N}, b::AbstractVector{N}) where {N<:Real}

Convert a matrix-vector representation to a linear-constraint representation.

Input

• A – matrix
• b – vector

Output

A list of linear constraints.

constrained_dimensions(hs::HalfSpace{N})::Vector{Int} where {N<:Real}

Return the indices in which a half-space is constrained.

Input

• hs – half-space

Output

A vector of ascending indices i such that the half-space is constrained in dimension i.

Examples

A 2D half-space with constraint $x1 ≥ 0$ is constrained in dimension 1 only.

translate(hs::HalfSpace{N}, v::AbstractVector{N}; share::Bool=false
         ) where {N<:Real}

Translate (i.e., shift) a half-space by a given vector.

Input

• hs – half-space
• v – translation vector
• share – (optional, default: false) flag for sharing unmodified parts of the original set representation

Output

A translated half-space.

Notes

The normal vectors of the halfspace (vector a in a⋅x ≤ b) is shared with the original halfspace if share == true.

Algorithm

A half-space $a⋅x ≤ b$ is transformed to the half-space $a⋅x ≤ b + a⋅v$. In other words, we add the dot product $a⋅v$ to $b$.

halfspace_left(p::AbstractVector{N},
               q::AbstractVector{N})::HalfSpace{N} where {N<:Real}

Return a half-space describing the 'left' of a two-dimensional line segment through two points.

Input

• p – first point
• q – second point

Output

The half-space whose boundary goes through the two points p and q and which describes the left-hand side of the directed line segment pq.

Algorithm

The implementation is simple: the half-space $a⋅x ≤ b$ is calculated as a = [dy, -dx], where $d = (dx, dy)$ denotes the line segment pq, that is, $\vec{d} = \vec{p} - \vec{q}$, and b = dot(p, a).

Examples

The left half-space of the "east" and "west" directions in two-dimensions are the upper and lower half-spaces:

julia> using LazySets: halfspace_left

julia> halfspace_left([0.0, 0.0], [1.0, 0.0])
HalfSpace{Float64,Array{Float64,1}}([0.0, -1.0], 0.0)

julia> halfspace_left([0.0, 0.0], [-1.0, 0.0])
HalfSpace{Float64,Array{Float64,1}}([0.0, 1.0], 0.0)

We create a box from the sequence of line segments that describe its edges:

julia> H1 = halfspace_left([-1.0, -1.0], [1.0, -1.0]);

julia> H2 = halfspace_left([1.0, -1.0], [1.0, 1.0]);

julia> H3 = halfspace_left([1.0, 1.0], [-1.0, 1.0]);

julia> H4 = halfspace_left([-1.0, 1.0], [-1.0, -1.0]);

julia> H = HPolygon([H1, H2, H3, H4]);

julia> B = BallInf([0.0, 0.0], 1.0);

julia> B ⊆ H && H ⊆ B
true

halfspace_right(p::AbstractVector{N},
                q::AbstractVector{N})::HalfSpace{N} where {N<:Real}

Return a half-space describing the 'right' of a two-dimensional line segment through two points.

Input

• p – first point
• q – second point

Output

The half-space whose boundary goes through the two points p and q and which describes the right-hand side of the directed line segment pq.

Algorithm

See the documentation of halfspace_left.

tosimplehrep(constraints::AbstractVector{LC})
    where {N<:Real, LC<:LinearConstraint{N}}

Return the simple H-representation $Ax ≤ b$ from a list of linear constraints.

Input

• constraints – a list of linear constraints

Output

The tuple (A, b) where A is the matrix of normal directions and b is the vector of offsets.

remove_redundant_constraints(constraints::AbstractVector{LC};
    backend=GLPKSolverLP()) where {N<:Real, LC<:LinearConstraint{N}}

Remove the redundant constraints of a given list of linear constraints.

Input

• constraints – list of constraints
• backend – (optional, default: GLPKSolverLP) numeric LP solver backend

Output

The list of constraints with the redundant ones removed, or an empty set if the constraints are infeasible.

Algorithm

See remove_redundant_constraints!(::AbstractVector{<:LinearConstraint{<:Real}}) for details.

 remove_redundant_constraints!(constraints::AbstractVector{LC};
     [backend]=GLPKSolverLP())::Bool where {N<:Real,
                                            LC<:LinearConstraint{N}}

Remove the redundant constraints of a given list of linear constraints; the list is updated in-place.

Input

• constraints – list of constraints
• backend – (optional, default: GLPKSolverLP) numeric LP solver backend

Output

true if the function was successful and the list of constraints constraints is modified by removing the redundant constraints, and false only if the constraints are infeasible.

Notes

Note that the result may be true even if the constraints are infeasible. For example, $x ≤ 0 && x ≥ 1$ will return true without removing any constraint. To check if the constraints are infeasible, use isempty(HPolyhedron(constraints).

Algorithm

If there are m constraints in n dimensions, this function checks one by one if each of the m constraints is implied by the remaining ones.

To check if the k-th constraint is redundant, an LP is formulated using the constraints that have not yet being removed. If, at an intermediate step, it is detected that a subgroup of the constraints is infeasible, this function returns false. If the calculation finished successfully, this function returns true.

For details, see Fukuda's Polyhedra FAQ.

Hyperplane{N<:Real} <: AbstractPolyhedron{N}

Type that represents a hyperplane of the form $a⋅x = b$.

Fields

• a – normal direction (non-zero)
• b – constraint

Examples

The plane $y = 0$:

julia> Hyperplane([0, 1.], 0.)
Hyperplane{Float64}([0.0, 1.0], 0.0)

dim(hp::Hyperplane)::Int

Return the dimension of a hyperplane.

Input

• hp – hyperplane

Output

The ambient dimension of the hyperplane.

ρ(d::AbstractVector{N}, hp::Hyperplane{N})::N where {N<:Real}

Evaluate the support function of a hyperplane in a given direction.

Input

• d – direction
• hp – hyperplane

Output

The support function of the hyperplane. If the set is unbounded in the given direction, the result is Inf.

σ(d::AbstractVector{N}, hp::Hyperplane{N}) where {N<:Real}

Return the support vector of a hyperplane.

Input

• d – direction
• hp – hyperplane

Output

The support vector in the given direction, which is only defined in the following two cases:

1. The direction has norm zero.
2. The direction is the hyperplane's normal direction or its opposite direction.

In all cases, the result is any point on the hyperplane. Otherwise this function throws an error.

∈(x::AbstractVector{N}, hp::Hyperplane{N})::Bool where {N<:Real}

Check whether a given point is contained in a hyperplane.

Input

• x – point/vector
• hp – hyperplane

Output

true iff $x ∈ hp$.

Algorithm

We just check if $x$ satisfies $a⋅x = b$.

an_element(hp::Hyperplane{N})::Vector{N} where {N<:Real}

Return some element of a hyperplane.

Input

• hp – hyperplane

Output

An element on the hyperplane.

rand(::Type{Hyperplane}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::Hyperplane{N}

Create a random hyperplane.

Input

• Hyperplane – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A random hyperplane.

Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1. Additionally, the constraint a is nonzero.

isbounded(hp::Hyperplane)::Bool

Determine whether a hyperplane is bounded.

Input

• hp – hyperplane

Output

false.

isuniversal(hp::Hyperplane{N}, [witness]::Bool=false
           )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether a hyperplane is universal.

Input

• P – hyperplane
• witness – (optional, default: false) compute a witness if activated

Output

• If witness option is deactivated: false
• If witness option is activated: (false, v) where $v ∉ P$

Algorithm

A witness is produced by adding the normal vector to an element on the hyperplane.

isempty(hp::Hyperplane)::Bool

Return if a hyperplane is empty or not.

Input

• hp – hyperplane

Output

false.

constrained_dimensions(hp::Hyperplane{N})::Vector{Int} where {N<:Real}

Return the indices in which a hyperplane is constrained.

Input

• hp – hyperplane

Output

A vector of ascending indices i such that the hyperplane is constrained in dimension i.

Examples

A 2D hyperplane with constraint $x1 = 0$ is constrained in dimension 1 only.

constraints_list(hp::Hyperplane{N}) where {N<:Real}

Return the list of constraints of a hyperplane.

Input

• hp – hyperplane

Output

A list containing two half-spaces.

translate(hp::Hyperplane{N}, v::AbstractVector{N}; share::Bool=false
         ) where {N<:Real}

Translate (i.e., shift) a hyperplane by a given vector.

Input

• hp – hyperplane
• v – translation vector
• share – (optional, default: false) flag for sharing unmodified parts of the original set representation

Output

A translated hyperplane.

Notes

The normal vectors of the hyperplane (vector a in a⋅x = b) is shared with the original hyperplane if share == true.

Algorithm

A hyperplane $a⋅x = b$ is transformed to the hyperplane $a⋅x = b + a⋅v$. In other words, we add the dot product $a⋅v$ to $b$.

Hyperrectangle{N<:Real} <: AbstractHyperrectangle{N}

Type that represents a hyperrectangle.

A hyperrectangle is the Cartesian product of one-dimensional intervals.

Fields

• center – center of the hyperrectangle as a real vector
• radius – radius of the ball as a real vector, i.e., half of its width along each coordinate direction

Examples

There is also a constructor from lower and upper bounds with keyword arguments high and low. The following two constructions are equivalent:

julia> c = ones(2);

julia> r = [0.1, 0.2];

julia> l = [0.9, 0.8];

julia> h = [1.1, 1.2];

julia> Hyperrectangle(c, r)
Hyperrectangle{Float64}([1.0, 1.0], [0.1, 0.2])
julia> Hyperrectangle(low=l, high=h)
Hyperrectangle{Float64}([1.0, 1.0], [0.1, 0.2])

rand(::Type{Hyperrectangle}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::Hyperrectangle{N}

Create a random hyperrectangle.

Input

• Hyperrectangle – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A random hyperrectangle.

Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1. Additionally, the radius is nonnegative.

center(H::Hyperrectangle{N})::Vector{N} where {N<:Real}

Return the center of a hyperrectangle.

Input

• H – hyperrectangle

Output

The center of the hyperrectangle.

radius_hyperrectangle(H::Hyperrectangle{N})::Vector{N} where {N<:Real}

Return the box radius of a hyperrectangle in every dimension.

Input

• H – hyperrectangle

Output

The box radius of the hyperrectangle.

radius_hyperrectangle(H::Hyperrectangle{N}, i::Int)::N where {N<:Real}

Return the box radius of a hyperrectangle in a given dimension.

Input

• H – hyperrectangle
• i – dimension of interest

Output

The radius in the given dimension.

translate(H::Hyperrectangle{N}, v::AbstractVector{N}; share::Bool=false
         ) where {N<:Real}

Translate (i.e., shift) a hyperrectangle by a given vector.

Input

• H – hyperrectangle
• v – translation vector
• share – (optional, default: false) flag for sharing unmodified parts of the original set representation

Output

A translated hyperrectangle.

Notes

The radius vector is shared with the original hyperrectangle if share == true.

Algorithm

We add the vector to the center of the hyperrectangle.

Interval{N<:Real, IN<:AbstractInterval{N}} <: AbstractHyperrectangle{N}

Type representing an interval on the real line. Mathematically, it is of the form

Fields

• dat – data container for the given interval

Notes

This type relies on the IntervalArithmetic.jl library for representation of intervals and arithmetic operations.

Examples

Unidimensional intervals are symbolic representations of a real closed interval.

We can create intervals in different ways, the simpler way is to pass a pair of numbers:

julia> x = Interval(0.0, 1.0)
Interval{Float64,IntervalArithmetic.Interval{Float64}}([0, 1])

or a 2-vector:

julia> x = Interval([0.0, 1.0])
Interval{Float64,IntervalArithmetic.Interval{Float64}}([0, 1])

Note that if the package IntervalArithmetic is loaded in the current scope, you have to prepend LazySets. to the Interval type since there is a name conflict otherwise.

julia> using IntervalArithmetic
WARNING: using IntervalArithmetic.Interval in module Main conflicts with an existing identifier.

julia> x = Interval(IntervalArithmetic.Interval(0.0, 1.0))
Interval{Float64,IntervalArithmetic.Interval{Float64}}([0, 1])

julia> dim(x)
1

julia> center(x)
1-element Array{Float64,1}:
 0.5

This type is such that the usual pairwise arithmetic operators +, -, * trigger the corresponding interval arithmetic backend method, and return a new Interval object. For the symbolic Minkowksi sum, use MinkowskiSum or ⊕.

Interval of other numeric types can be created as well, eg. a rational interval:

julia> Interval(0//1, 2//1)
Interval{Rational{Int64},AbstractInterval{Rational{Int64}}}([0//1, 2//1])

dim(x::Interval)::Int

Return the ambient dimension of an interval.

Input

• x – interval

Output

The integer 1.

σ(d::AbstractVector{N}, x::Interval{N}) where {N<:Real}

Return the support vector of an interval in a given direction.

Input

• d – direction
• x – interval

Output

Support vector in the given direction.

ρ(d::AbstractVector{N}, x::Interval{N}) where {N<:Real}

Evaluate the support function of an interval in a given direction.

Input

• d – direction
• x – interval

Output

Evaluation of the support function in the given direction.

∈(v::AbstractVector{N}, x::Interval{N}) where {N<:Real})

Return whether a vector is contained in the interval.

Input

• v – one-dimensional vector
• x – interval

Output

true iff x contains v's first component.

∈(v::N, x::Interval{N}) where {N<:Real}

Return whether a number is contained in the interval.

Input

• v – scalar
• x – interval

Output

true iff x contains v.

an_element(x::Interval{N})::Vector{N} where {N<:Real}

Return some element of an interval.

Input

• x – interval

Output

The left border (min(x)) of the interval.

vertices_list(x::Interval{N})::Vector{Vector{N}} where {N<:Real}

Return the list of vertices of this interval.

Input

• x – interval

Output

The list of vertices of the interval represented as two one-dimensional vectors.

translate(x::Interval{N}, v::AbstractVector{N}) where {N<:Real}

Translate (i.e., shift) an interval by a given vector.

Input

• x – interval
• v – translation vector

Output

A translated interval.

Algorithm

We add the vector to the left and right of the interval.

center(x::Interval{N})::Vector{N} where {N<:Real}

Return the interval's center.

Input

• x – interval

Output

The center, or midpoint, of $x$.

min(x::Interval{N})::N where {N<:Real}

Return the lower component of an interval.

Input

• x – interval

Output

The lower (lo) component of the interval.

max(x::Interval{N})::N where {N<:Real}

Return the higher or upper component of an interval.

Input

• x – interval

Output

The higher (hi) component of the interval.

low(x::Interval{N})::Vector{N} where {N<:Real}

Return the lower coordinate of an interval set.

Input

• x – interval

Output

A vector with the lower coordinate of the interval.

high(x::Interval{N})::Vector{N} where {N<:Real}

Return the higher coordinate of an interval set.

Input

• x – interval

Output

A vector with the higher coordinate of the interval.

radius_hyperrectangle(x::Interval{N})::Vector{N} where {N<:Real}

Return the box radius of an interval in every dimension.

Input

• x – interval

Output

The box radius of the interval (a one-dimensional vector).

radius_hyperrectangle(x::Interval{N}, i::Int)::N where {N<:Real}

Return the box radius of an interval in a given dimension.

Input

• x – interval
• i – dimension index (must be 1)

Output

The box radius in the given dimension.

+(x::Interval{N}, y::Interval{N}) where {N<:Real}

Return the sum of the intervals.

Input

• x – interval
• y – interval

Output

The sum of the intervals as a new Interval set.

-(x::Interval{N}, y::Interval{N}) where {N<:Real}

Return the difference of the intervals.

Input

• x – interval
• y – interval

Output

The difference of the intervals as a new Interval set.

    *(x::Interval{N}, y::Interval{N}) where {N<:Real}

Return the product of the intervals.

Input

• x – interval
• y – interval

Output

The product of the intervals as a new Interval set.

rand(::Type{Interval}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::Interval{N}

Create a random interval.

Input

• Interval – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 1) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A random interval.

Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1.

isflat(I::Interval)::Bool

Determine whether an interval is flat, i.e. whether its extreme values coincide.

Input

• I – interval

Output

A boolean which is true if the interval is flat and false otherwise.

Notes

For robustness with respect to floating-point inputs, this function relies on the result of isapproxzero when applied to the diameter of the interval. Hence, this function depends on the absolute zero tolerance ABSZTOL.

plot_recipe(I::Interval{N}, [ε]::N=zero(N)) where {N<:Real}

Convert an interval to a pair (x, y) of points for plotting.

Input

• I – interval
• ε – (optional, default: 0) ignored, used for dispatch

Output

A pair (x, y) of two points that can be plotted.

Notes

We consider the interval as a line segment with y coordinate equal to zero.

Line{N<:Real, VN<:AbstractVector{N}} <: AbstractPolyhedron{N}

Type that represents a line in 2D of the form $a⋅x = b$ (i.e., a special case of a Hyperplane).

Fields

• a – normal direction (non-zero)
• b – constraint

Examples

The line $y = -x + 1$:

julia> Line([1., 1.], 1.)
Line{Float64,Array{Float64,1}}([1.0, 1.0], 1.0)

dim(L::Line)::Int

Return the ambient dimension of a line.

Input

• L – line

Output

The ambient dimension of the line, which is 2.

σ(d::AbstractVector{N}, L::Line{N}) where {N<:Real}

Return the support vector of a line in a given direction.

Input

• d – direction
• L – line

Output

The support vector in the given direction, which is defined the same way as for the more general Hyperplane.

∈(x::AbstractVector{N}, L::Line{N})::Bool where {N<:Real}

Check whether a given point is contained in a line.

Input

• x – point/vector
• L – line

Output

true iff x ∈ L.

Algorithm

The point $x$ belongs to the line if and only if $a⋅x = b$ holds.

an_element(L::Line{N})::Vector{N} where {N<:Real}

Return some element of a line.

Input

• L – line

Output

An element on the line.

Algorithm

If the $b$ value of the line is zero, the result is the origin. Otherwise the result is some $x = [x1, x2]$ such that $a·[x1, x2] = b$. We first find out in which dimension $a$ is nonzero, say, dimension 1, and then choose $x1 = 1$ and accordingly $x2 = \frac{b - a1}{a2}$.

rand(::Type{Line}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::Line{N}

Create a random line.

Input

• Line – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A random line.

Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1. Additionally, the constraint a is nonzero.

isbounded(L::Line)::Bool

Determine whether a line is bounded.

Input

• L – line

Output

false.

isuniversal(L::Line{N}, [witness]::Bool=false
           )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether a line is universal.

Input

• P – line
• witness – (optional, default: false) compute a witness if activated

Output

• If witness option is deactivated: false
• If witness option is activated: (false, v) where $v ∉ P$

Algorithm

Witness production falls back to isuniversal(::Hyperplane).

isempty(L::Line)::Bool

Return if a line is empty or not.

Input

• L – line

Output

false.

constrained_dimensions(L::Line{N})::Vector{Int} where {N<:Real}

Return the indices in which a line is constrained.

Input

• L – line

Output

A vector of ascending indices i such that the line is constrained in dimension i.

Examples

A line with constraint $x1 = 0$ is constrained in dimension 1 only.

constraints_list(L::Line{N}) where {N<:Real}

Return the list of constraints of a line.

Input

• L – line

Output

A list containing two half-spaces.

translate(L::Line{N}, v::AbstractVector{N}; share::Bool=false
         ) where {N<:Real}

Translate (i.e., shift) a line by a given vector.

Input

• L – line
• v – translation vector
• share – (optional, default: false) flag for sharing unmodified parts of the original set representation

Output

A translated line.

Notes

The normal vector of the line (vector a in a⋅x = b) is shared with the original line if share == true.

Algorithm

A line $a⋅x = b$ is transformed to the line $a⋅x = b + a⋅v$. In other words, we add the dot product $a⋅v$ to $b$.

LineSegment{N<:Real} <: AbstractZonotope{N}

Type that represents a line segment in 2D between two points $p$ and $q$.

Fields

• p – first point
• q – second point

Examples

A line segment along the $x = y$ diagonal:

julia> s = LineSegment([0., 0], [1., 1.])
LineSegment{Float64}([0.0, 0.0], [1.0, 1.0])
julia> dim(s)
2

Use plot(s) to plot the extreme points of s and the line segment joining them. Membership test is computed with ∈ (in):

julia> [0., 0] ∈ s && [.25, .25] ∈ s && [1., 1] ∈ s && [.5, .25] ∉ s
true

We can check the intersection with another line segment, and optionally compute a witness (which is just the common point in this case):

julia> sn = LineSegment([1., 0], [0., 1.])
LineSegment{Float64}([1.0, 0.0], [0.0, 1.0])
julia> isempty(s ∩ sn)
false
julia> is_intersection_empty(s, sn, true)
(false, [0.5, 0.5])

dim(L::LineSegment)::Int

Return the ambient dimension of a line segment.

Input

• L – line segment

Output

The ambient dimension of the line segment, which is 2.

σ(d::AbstractVector{N}, L::LineSegment{N}) where {N<:Real}

Return the support vector of a line segment in a given direction.

Input

• d – direction
• L – line segment

Output

The support vector in the given direction.

Algorithm

If the angle between the vector $q - p$ and $d$ is bigger than 90° and less than 270° (measured in counter-clockwise order), the result is $p$, otherwise it is $q$. If the angle is exactly 90° or 270°, or if the direction has norm zero, this implementation returns $q$.

∈(x::AbstractVector{N}, L::LineSegment{N})::Bool where {N<:Real}

Check whether a given point is contained in a line segment.

Input

• x – point/vector
• L – line segment

Output

true iff $x ∈ L$.

Algorithm

Let $L = (p, q)$ be the line segment with extremes $p$ and $q$, and let $x$ be the given point.

1. A necessary condition for $x ∈ (p, q)$ is that the three points are aligned, thus their cross product should be zero.
2. It remains to check that $x$ belongs to the box approximation of $L$. This amounts to comparing each coordinate with those of the extremes $p$ and $q$.

Notes

The algorithm is inspired from here.

center(L::LineSegment{N})::Vector{N} where {N<:Real}

Return the center of a line segment.

Input

• L – line segment

Output

The center of the line segment.

an_element(L::LineSegment{N}) where {N<:Real}

Return some element of a line segment.

Input

• L – line segment

Output

The first vertex of the line segment.

rand(::Type{LineSegment}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::LineSegment{N}

Create a random line segment.

Input

• LineSegment – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A random line segment.

Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1.

halfspace_left(L::LineSegment)

Return a half-space describing the 'left' of a two-dimensional line segment through two points.

Input

• L – line segment

Output

The half-space whose boundary goes through the two points p and q and which describes the left-hand side of the directed line segment pq.

halfspace_right(L::LineSegment)

Return a half-space describing the 'right' of a two-dimensional line segment through two points.

Input

• L – line segment

Output

The half-space whose boundary goes through the two points p and q and which describes the right-hand side of the directed line segment pq.

vertices_list(L::LineSegment{N}
             )::Vector{<:AbstractVector{N}} where {N<:Real}

Return the list of vertices of a line segment.

Input

• L – line segment

Output

The list of end points of the line segment.

constraints_list(L::LineSegment{N}) where {N<:Real}

Return the list of constraints defining a line segment in 2D.

Input

• L – line segment

Output

A vector of constraints that define the line segment.

Algorithm

$L$ is defined by 4 constraints. In this algorithm, the first two constraints are returned by $halfspace_right$ and $halfspace_left$, and the other two are obtained by considering the vector normal to the line segment that passes through each opposite vertex.

Notes

This function returns a vector of halfspaces. It does not return equality constraints.

translate(L::LineSegment{N}, v::AbstractVector{N}) where {N<:Real}

Translate (i.e., shift) a line segment by a given vector.

Input

• L – line segment
• v – translation vector

Output

A translated line segment.

Algorithm

We add the vector to both defining points of the line segment.

generators(L::LineSegment{N}) where {N<:Real}

Return an iterator over the (single) generator of a line segment.

Input

• L – line segment

Output

A one-element iterator over the generator of L.

genmat(L::LineSegment)

Return the generator matrix of a line segment.

Input

• L – line segment

Output

A matrix with a single column representing the generator of L.

plot_recipe(L::LineSegment{N}, [ε]::N=zero(N)) where {N<:Real}

Convert a line segment to a pair (x, y) of points for plotting.

Input

• L – line segment
• ε – (optional, default: 0) ignored, used for dispatch

Output

A pair (x, y) of two points that can be plotted.

plot_linesegment(X::Union{Interval{N}, LineSegment{N}}, [ε]::N=zero(N); ...)
    where {N<:Real}

Plot a line segment or an interval.

Input

• X – line segment or interval
• ε – (optional, default: 0) ignored, used for dispatch

Examples

julia> L = LineSegment([0., 0.], [1., 1.]);

julia> plot(L)

To control the color of the line, use the linecolor keyword argument, and to control the color of the end points, use the markercolor keyword argument. To control the width, use linewidth.

julia> plot(L, markercolor="green", linecolor="red", linewidth=2.)

To omit the markers, use markershape=:none. You also need to pass a value for seriestype=:path explicitly (this seems to be an external bug).

julia> plot(L, seriestype=:path, markershape=:none)

A shorter alternative is to pass marker=0, but this may result in small dots as markers based on the plotting backend.

julia> plot(L, marker=0)

HPolygon{N<:Real} <: AbstractHPolygon{N}

Type that represents a convex polygon in constraint representation whose edges are sorted in counter-clockwise fashion with respect to their normal directions.

Fields

• constraints – list of linear constraints, sorted by the normal direction in counter-clockwise fashion
• sort_constraints – (optional, default: true) flag for sorting the constraints (sortedness is a running assumption of this type)
• check_boundedness – (optional, default: false) flag for checking if the constraints make the polygon bounded; (boundedness is a running assumption of this type)
• prune – (optional, default: true) flag for removing redundant constraints

Notes

The option sort_constraints can be used to deactivate automatic sorting of constraints in counter-clockwise fashion, which is an invariant of this type. Alternatively, one can construct an HPolygon with empty constraints list, which can then be filled iteratively using addconstraint!.

Similarly, the option prune can be used to deactivate automatic pruning of redundant constraints.

Another type assumption is that the polygon is bounded. The option check_boundedness can be used to assert this. This option is deactivated by default because we explicitly want to allow the iterative addition of the constraints, and hence one has to initially construct an empty list of constraints (which represents an unbounded set). The user has to make sure that the HPolygon is not used before the constraints actually describe a bounded set. The function isbounded can be used to manually assert boundedness.

σ(d::AbstractVector{N}, P::HPolygon{N};
  [linear_search]::Bool=(length(P.constraints) < BINARY_SEARCH_THRESHOLD)
 ) where {N<:Real}

Return the support vector of a polygon in a given direction.

Input

• d – direction
• P – polygon in constraint representation
• linear_search – (optional, default: see below) flag for controlling whether to perform a linear search or a binary search

Output

The support vector in the given direction. The result is always one of the vertices; in particular, if the direction has norm zero, any vertex is returned.

Algorithm

Comparison of directions is performed using polar angles; see the overload of <= for two-dimensional vectors.

For polygons with BINARY_SEARCH_THRESHOLD = 10 or more constraints we use a binary search by default.

translate(v::AbstractVector{N}, P::HPolygon{N}; share::Bool=false
         ) where {N<:Real}

Translate (i.e., shift) a polygon in constraint representation by a given vector.

Input

• P – polygon in constraint representation
• v – translation vector
• share – (optional, default: false) flag for sharing unmodified parts of the original set representation

Output

A translated polygon in constraint representation.

Notes

The normal vectors of the constraints (vector a in a⋅x ≤ b) are shared with the original constraints if share == true.

Algorithm

We translate every constraint.

HPolygonOpt{N<:Real} <: AbstractHPolygon{N}

Type that represents a convex polygon in constraint representation whose edges are sorted in counter-clockwise fashion with respect to their normal directions. This is a refined version of HPolygon.

Fields

• constraints – list of linear constraints, sorted by the normal direction in counter-clockwise fashion
• ind – index in the list of constraints to begin the search to evaluate the support function
• sort_constraints – (optional, default: true) flag for sorting the constraints (sortedness is a running assumption of this type)
• check_boundedness – (optional, default: false) flag for checking if the constraints make the polygon bounded; (boundedness is a running assumption of this type)
• prune – (optional, default: true) flag for removing redundant constraints

Notes

This structure is optimized to evaluate the support function/vector with a large sequence of directions that are close to each other. The strategy is to have an index that can be used to warm-start the search for optimal values in the support vector computation.

The option sort_constraints can be used to deactivate automatic sorting of constraints in counter-clockwise fashion, which is an invariant of this type. Alternatively, one can construct an HPolygonOpt with empty constraints list, which can then be filled iteratively using addconstraint!.

Similarly, the option prune can be used to deactivate automatic pruning of redundant constraints.

Another type assumption is that the polygon is bounded. The option check_boundedness can be used to assert this. This option is deactivated by default because we explicitly want to allow the iterative addition of the constraints, and hence one has to initially construct an empty list of constraints (which represents an unbounded set). The user has to make sure that the HPolygonOpt is not used before the constraints actually describe a bounded set. The function isbounded can be used to manually assert boundedness.

σ(d::AbstractVector{N}, P::HPolygonOpt{N};
  [linear_search]::Bool=(length(P.constraints) < BINARY_SEARCH_THRESHOLD)
 ) where {N<:Real}

Return the support vector of an optimized polygon in a given direction.

Input

• d – direction
• P – optimized polygon in constraint representation
• linear_search – (optional, default: see below) flag for controlling whether to perform a linear search or a binary search

Output

The support vector in the given direction. The result is always one of the vertices; in particular, if the direction has norm zero, any vertex is returned.

Algorithm

Comparison of directions is performed using polar angles; see the overload of <= for two-dimensional vectors.

For polygons with BINARY_SEARCH_THRESHOLD = 10 or more constraints we use a binary search by default.

translate(P::HPolygonOpt{N}, v::AbstractVector{N}; share::Bool=false
         ) where {N<:Real}

Translate (i.e., shift) an optimized polygon in constraint representation by a given vector.

Input

• P – optimized polygon in constraint representation
• v – translation vector
• share – (optional, default: false) flag for sharing unmodified parts of the original set representation

Output

A translated optimized polygon in constraint representation.

Notes

The normal vectors of the constraints (vector a in a⋅x ≤ b) are shared with the original constraints if share == true.

Algorithm

We translate every constraint.

VPolygon{N<:Real, VN<:AbstractVector{N}} <: AbstractPolygon{N}

Type that represents a polygon by its vertices.

Fields

• vertices – the list of vertices

Notes

The constructor of VPolygon runs a convex hull algorithm on its vertices by default, to remove the possibly redundant vertices. The vertices are sorted in counter-clockwise fashion. Use the flag apply_convex_hull=false to skip the computation of the convex hull.

• VPolygon(vertices::Vector{Vector{N}}; apply_convex_hull::Bool=true, algorithm::String="monotone_chain")

σ(d::AbstractVector{N}, P::VPolygon{N}) where {N<:Real}

Return the support vector of a polygon in a given direction.

Input

• d – direction
• P – polygon in vertex representation

Output

The support vector in the given direction. If the direction has norm zero, the first vertex is returned.

Algorithm

This implementation uses a binary search algorithm when the polygon has more than ten vertices and a brute-force search when it has ten or less. For the brute-force search, it compares the projection of each vector along the given direction and runs in $O(n)$ where $n$ is the number of vertices. For the binary search the algorithm runs in $O(log n)$. We follow this implementation based on an algorithm described in [1].

[1] Joseph O'Rourke, Computational Geometry in C (2nd Edition)

∈(x::AbstractVector{N}, P::VPolygon{N})::Bool where {N<:Real}

Check whether a given point is contained in a polygon in vertex representation.

Input

• x – point/vector
• P – polygon in vertex representation

Output

true iff $x ∈ P$.

Algorithm

This implementation exploits that the polygon's vertices are sorted in counter-clockwise fashion. Under this assumption we can just check if the vertex lies on the left of each edge, using the dot product.

Examples

julia> P = VPolygon([[2.0, 3.0], [3.0, 1.0], [5.0, 1.0], [4.0, 5.0]];
                    apply_convex_hull=false);

julia> [4.5, 3.1] ∈ P
false
julia> [4.5, 3.0] ∈ P
true
julia> [4.4, 3.4] ∈ P  #  point lies on the edge -> floating-point error
false
julia> P = VPolygon([[2//1, 3//1], [3//1, 1//1], [5//1, 1//1], [4//1, 5//1]];
                     apply_convex_hull=false);

julia> [44//10, 34//10] ∈ P  #  with rational numbers the answer is correct
true

an_element(P::VPolygon{N})::Vector{N} where {N<:Real}

Return some element of a polygon in vertex representation.

Input

• P – polygon in vertex representation

Output

The first vertex of the polygon in vertex representation.

rand(::Type{VPolygon}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::VPolygon{N}

Create a random polygon in vertex representation.

Input

• VPolygon – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding
• num_vertices – (optional, default: -1) number of vertices of the polygon (see comment below)

Output

A random polygon in vertex representation.

Algorithm

We follow the idea here based on P. Valtr. Probability that n random points are in convex position. There is also a nice video available here.

The number of vertices can be controlled with the argument num_vertices. For a negative value we choose a random number in the range 3:10.

vertices_list(P::VPolygon{N})::Vector{Vector{N}} where {N<:Real}

Return the list of vertices of a convex polygon in vertex representation.

Input

• P – a polygon vertex representation

Output

List of vertices.

tohrep(P::VPolygon{N}, ::Type{HPOLYGON}=HPolygon
      ) where {N<:Real, HPOLYGON<:AbstractHPolygon}

Build a constraint representation of the given polygon.

Input

• P – polygon in vertex representation
• HPOLYGON – (optional, default: HPolygon) type of target polygon

Output

The same polygon but in constraint representation, an AbstractHPolygon.

Algorithm

The algorithms consists of adding an edge for each consecutive pair of vertices. Since the vertices are already ordered in counter-clockwise fashion (CWW), the constraints will be sorted automatically (CCW) if we start with the first edge between the first and second vertex.

tovrep(P::VPolygon{N})::VPolygon{N} where {N<:Real}

Build a vertex representation of the given polygon.

Input

• P – polygon in vertex representation

Output

The identity, i.e., the same polygon instance.

constraints_list(P::VPolygon{N}) where {N<:Real}

Return the list of constraints defining a polygon in V-representation.

Input

• P – polygon in V-representation

Output

The list of constraints of the polygon.

Algorithm

First the H-representation of $P$ is computed, then its list of constraints is returned.

translate(P::VPolygon{N}, v::AbstractVector{N}) where {N<:Real}

Translate (i.e., shift) a polygon in vertex representation by a given vector.

Input

• P – polygon in vertex representation
• v – translation vector

Output

A translated polygon in vertex representation.

Algorithm

We add the vector to each vertex of the polygon.

remove_redundant_vertices(P::VPolygon{N};
                          [algorithm]::String="monotone_chain")::VPolygon{N} where {N<:Real}

Return the polygon obtained by removing the redundant vertices of the given polygon.

Input

• P – polygon in vertex representation
• algorithm – (optional, default: "monotone_chain") the algorithm used to compute the convex hull

Output

A new polygon such that its vertices are the convex hull of the given polygon.

Algorithm

A convex hull algorithm is used to compute the convex hull of the vertices of the given input polygon P; see ?convex_hull for details on the available algorithms. The vertices of the output polygon are sorted in counter-clockwise fashion.

remove_redundant_vertices!(P::VPolygon{N};
                           [algorithm]::String="monotone_chain")::VPolygon{N} where {N<:Real}

Remove the redundant vertices of the given polygon.

Input

• P – polygon in vertex representation
• algorithm – (optional, default: "monotone_chain") the algorithm used to compute the convex hull

Output

A new polygon such that its vertices are the convex hull of the given polygon.

Algorithm

A convex hull algorithm is used to compute the convex hull of the vertices of the given input polygon P; see ?convex_hull for details on the available algorithms. The vertices of the output polygon are sorted in counter-clockwise fashion.

jump2pi(x::N)::N where {N<:AbstractFloat}

Return $x + 2π$ if $x$ is negative, otherwise return $x$.

Input

• x – real scalar

Output

$x + 2π$ if $x$ is negative, $x$ otherwise.

Examples

julia> using LazySets: jump2pi

julia> jump2pi(0.0)
0.0

julia> jump2pi(-0.5)
5.783185307179586

julia> jump2pi(0.5)
0.5

<=(u::AbstractVector{N}, v::AbstractVector{N})::Bool where {N<:AbstractFloat}

Compares two 2D vectors by their direction.

Input

• u – first 2D direction
• v – second 2D direction

Output

True iff $\arg(u) [2π] ≤ \arg(v) [2π]$

Notes

The argument is measured in counter-clockwise fashion, with the 0 being the direction (1, 0).

Algorithm

The implementation uses the arctangent function with sign, atan, which for two arguments implements the atan2 function.

<=(u::AbstractVector{N}, v::AbstractVector{N})::Bool where {N<:Real}

Compare two 2D vectors by their direction.

Input

• u – first 2D direction
• v – second 2D direction

Output

true iff $\arg(u) [2π] ≤ \arg(v) [2π]$.

Notes

The argument is measured in counter-clockwise fashion, with the 0 being the direction (1, 0).

Algorithm

The implementation checks the quadrant of each direction, and compares directions using the right-hand rule (see is_right_turn). In particular, this method does not use the arctangent.

quadrant(w::AbstractVector{N})::Int where {N<:Real}

Compute the quadrant where the direction w belongs.

Input

• w – direction

Output

An integer from 0 to 3, with the following convention:

     ^
   1 | 0
  ---+-->
   2 | 3

Algorithm

The idea is to encode the following logic function: $11 ↦ 0, 01 ↦ 1, 00 ↦ 2, 10 ↦ 3$, according to the convention of above.

This function is inspired from AGPX's answer in: Sort points in clockwise order?

HPolytope{N<:Real} <: AbstractPolytope{N}

Type that represents a convex polytope in H-representation.

Fields

• constraints – vector of linear constraints
• check_boundedness – (optional, default: false) flag for checking if the constraints make the polytope bounded; (boundedness is a running assumption of this type)

Note

Recall that a polytope is a bounded polyhedron. Boundedness is a running assumption in this type.

HPolyhedron{N<:Real} <: AbstractPolyhedron{N}

Type that represents a convex polyhedron in H-representation.

Fields

• constraints – vector of linear constraints

dim(P::HPoly{N})::Int where {N<:Real}

Return the dimension of a polyhedron in H-representation.

Input

• P – polyhedron in H-representation

Output

The ambient dimension of the polyhedron in H-representation. If it has no constraints, the result is $-1$.

ρ(d::AbstractVector{N}, P::HPoly{N})::N where {N<:Real}

Evaluate the support function of a polyhedron (in H-representation) in a given direction.

Input

• d – direction
• P – polyhedron in H-representation

Output

The support function of the polyhedron. If a polytope is unbounded in the given direction, we throw an error. If a polyhedron is unbounded in the given direction, the result is Inf.

Algorithm

This implementation uses GLPKSolverLP as linear programming backend.

σ(d::AbstractVector{N}, P::HPoly{N}) where {N<:Real}

Return the support vector of a polyhedron (in H-representation) in a given direction.

Input

• d – direction
• P – polyhedron in H-representation

Output

The support vector in the given direction.

Algorithm

This implementation uses GLPKSolverLP as linear programming backend.

addconstraint!(P::HPoly{N},
               constraint::LinearConstraint{N})::Nothing where {N<:Real}

Add a linear constraint to a polyhedron in H-representation.

Input

• P – polyhedron in H-representation
• constraint – linear constraint to add

Output

Nothing.

Notes

It is left to the user to guarantee that the dimension of all linear constraints is the same.

constraints_list(P::HPoly{N}) where {N<:Real}

Return the list of constraints defining a polyhedron in H-representation.

Input

• P – polyhedron in H-representation

Output

The list of constraints of the polyhedron.

tohrep(P::HPoly{N}) where {N<:Real}

Return a constraint representation of the given polyhedron in constraint representation (no-op).

Input

• P – polyhedron in constraint representation

Output

The same polyhedron instance.

tovrep(P::HPoly{N};
      [backend]=default_polyhedra_backend(P, N)) where {N<:Real}

Transform a polyhedron in H-representation to a polytope in V-representation.

Input

• P – polyhedron in constraint representation
• backend – (optional, default: default_polyhedra_backend(P, N)) the backend for polyhedral computations

Output

The VPolytope which is the vertex representation of the given polyhedron in constraint representation.

Notes

The conversion may not preserve the numeric type (e.g., with N == Float32) depending on the backend. For further information on the supported backends see Polyhedra's documentation.

isempty(P::HPoly{N}, witness::Bool=false; [usepolyhedrainterface]::Bool=false, [solver]=GLPKSolverLP(), [backend]=nothing )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Determine whether a polyhedron is empty.

Input

• P – polyhedron
• witness – (optional, default: false) compute a witness if activated
• use_polyhedra_interface – (optional, default: false) if true, we use the Polyhedra interface for the emptiness test
• solver – (optional, default: GLPKSolverLP()) LP-solver backend
• backend – (optional, default: nothing) backend for polyhedral computations in Polyhedra; its value is set internally (see the Notes below for details)

Output

• If witness option is deactivated: true iff $P = ∅$
• If witness option is activated:
  • (true, []) iff $P = ∅$
  • (false, v) iff $P ≠ ∅$ and $v ∈ P$

Notes

The default value of the backend is set internally and depends on whether the use_polyhedra_interface option is set or not. If the option is set, we use default_polyhedra_backend(P, N).

Witness production is not supported if use_polyhedra_interface is true.

Algorithm

The algorithm sets up a feasibility LP for the constraints of P. If use_polyhedra_interface is true, we call Polyhedra.isempty. Otherwise, we set up the LP internally.

translate(P::PT, v::AbstractVector{N}; share::Bool=false
         ) where {N<:Real, PT<:HPoly{N}}

Translate (i.e., shift) a polyhedron in constraint representation by a given vector.

Input

• P – polyhedron in constraint representation
• v – translation vector
• share – (optional, default: false) flag for sharing unmodified parts of the original set representation

Output

A translated polyhedron in constraint representation.

Notes

The normal vectors of the constraints (vector a in a⋅x ≤ b) are shared with the original constraints if share == true.

Algorithm

We translate every constraint.

polyhedron(P::HPoly{N};
           [backend]=default_polyhedra_backend(P, N)) where {N<:Real}

Return an HRep polyhedron from Polyhedra.jl given a polytope in H-representation.

Input

• P – polytope
• backend – (optional, default: call default_polyhedra_backend(P, N)) the polyhedral computations backend

Output

An HRep polyhedron.

Notes

For further information on the supported backends see Polyhedra's documentation.

remove_redundant_constraints(P::PT;
                             backend=GLPKSolverLP()
                            )::Union{PT, EmptySet{N}} where {N<:Real,
                                                             PT<:HPoly{N}}

Remove the redundant constraints in a polyhedron in H-representation.

Input

• P – polyhedron
• backend – (optional, default: GLPKSolverLP) the numeric LP solver backend

Output

A polyhedron equivalent to P but with no redundant constraints, or an empty set if P is detected to be empty, which may happen if the constraints are infeasible.

Algorithm

See remove_redundant_constraints!(::Vector{LinearConstraint{<:Real}}) for details.

remove_redundant_constraints!(P::HPoly{N};
                              backend=GLPKSolverLP())::Bool where {N<:Real}

Remove the redundant constraints in a polyhedron in H-representation; the polyhedron is updated in-place.

Input

• P – polyhedron
• backend – (optional, default: GLPKSolverLP) the numeric LP solver backend

Output

true if the method was successful and the polyhedron P is modified by removing its redundant constraints, and false if P is detected to be empty, which may happen if the constraints are infeasible.

Algorithm

See remove_redundant_constraints!(::Vector{LinearConstraint{<:Real}}) for details.

rand(::Type{HPolytope}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::HPolytope{N}

Create a random polytope in constraint representation.

Input

• HPolytope – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding
• num_vertices – (optional, default: -1) upper bound on the number of vertices of the polytope (see comment below)

Output

A random polytope in constraint representation.

Algorithm

We create a random polytope in vertex representation and convert it to constraint representation (hence the argument num_vertices). See rand(::Type{VPolytope}).

vertices_list(P::HPolytope{N};
              [backend]=nothing,
              [prune]::Bool=true)::Vector{Vector{N}} where {N<:Real}

Return the list of vertices of a polytope in constraint representation.

Input

• P – polytope in constraint representation
• backend – (optional, default: nothing) the polyhedral computations backend
• prune – (optional, default: true) flag to remove redundant vertices

Output

List of vertices.

Algorithm

If the polytope is two-dimensional, the polytope is converted to a polygon in H-representation and then its vertices_list function is used. This ensures that, by default, the optimized two-dimensional methods are used.

It is possible to use the Polyhedra backend in two-dimensions as well by passing, e.g. backend=CDDLib.Library().

If the polytope is not two-dimensional, the concrete polyhedra manipulation library Polyhedra is used. The actual computation is performed by a given backend; for the default backend used in LazySets see default_polyhedra_backend(N). For further information on the supported backends see Polyhedra's documentation.

isbounded(P::HPolytope, [use_type_assumption]::Bool=true)::Bool

Determine whether a polytope in constraint representation is bounded.

Input

• P – polytope in constraint representation
• use_type_assumption – (optional, default: true) flag for ignoring the type assumption that polytopes are bounded

Output

true if use_type_assumption is activated. Otherwise, true iff P is bounded.

Algorithm

If !use_type_assumption, we convert P to an HPolyhedron P2 and then use isbounded(P2).

rand(::Type{HPolyhedron}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::HPolyhedron{N}

Create a polyhedron.

Input

• HPolyhedron – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension (is ignored)
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A polyhedron.

Algorithm

We first create a random polytope and then randomly remove some of the constraints.

isbounded(P::HPolyhedron)::Bool

Determine whether a polyhedron in constraint representation is bounded.

Input

• P – polyhedron in constraint representation

Output

true iff the polyhedron is bounded.

Algorithm

We first check if the polyhedron has more than max(dim(P), 1) constraints, which is a necessary condition for boundedness. If so, we check boundedness via isbounded_unit_dimensions.

isuniversal(P::HPolyhedron{N}, [witness]::Bool=false
           )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether a polyhedron is universal.

Input

• P – polyhedron
• witness – (optional, default: false) compute a witness if activated

Output

• If witness option is deactivated: true iff $P$ is universal
• If witness option is activated:
  • (true, []) iff $P$ is universal
  • (false, v) iff $P$ is not universal and $v ∉ P$

Algorithm

P is universal iff it has no constraints.

A witness is produced using isuniversal(H) where H is the first linear constraint of P.

vertices_list(P::HPolyhedron{N}) where {N<:Real}

Return the list of vertices of a polyhedron in constraint representation.

Input

• P – polyhedron in constraint representation

Output

This function returns an error because the polyhedron is possibly unbounded. If P is known to be bounded, try converting to HPolytope first:

julia> P = HPolyhedron([HalfSpace([1.0, 0.0], 1.0),
                        HalfSpace([0.0, 1.0], 1.0),
                        HalfSpace([-1.0, 0.0], 1.0),
                        HalfSpace([0.0, -1.0], 1.0)]);

julia> P_as_polytope = convert(HPolytope, P);

singleton_list(P::HPolyhedron{N})::Vector{Singleton{N}} where {N<:Real}

Return the vertices of a polyhedron in H-representation as a list of singletons.

Input

• P – polytope in constraint representation

Output

This function returns an error because the polyhedron is possibly unbounded. If P is known to be bounded, try converting to HPolytope first.

VPolytope{N<:Real} <: AbstractPolytope{N}

Type that represents a convex polytope in V-representation.

Fields

• vertices – the list of vertices

dim(P::VPolytope)::Int

Return the dimension of a polytope in V-representation.

Input

• P – polytope in V-representation

Output

The ambient dimension of the polytope in V-representation. If it is empty, the result is $-1$.

Examples

julia> v = VPolytope();

julia> dim(v) > 0
false

julia> v = VPolytope([ones(3)])
VPolytope{Float64}(Array{Float64,1}[[1.0, 1.0, 1.0]])

julia> dim(v) == 3
true

σ(d::AbstractVector{N}, P::VPolytope{N}) where {N<:Real}

Return the support vector of a polytope in V-representation in a given direction.

Input

• d – direction
• P – polytope in V-representation

Output

The support vector in the given direction.

Algorithm

A support vector maximizes the support function. For a polytope, the support function is always maximized in some vertex. Hence it is sufficient to check all vertices.

∈(x::AbstractVector{N}, P::VPolytope{N};
  solver=GLPKSolverLP(method=:Simplex))::Bool where {N<:Real}

Check whether a given point is contained in a polytope in vertex representation.

Input

• x – point/vector
• P – polytope in vertex representation

Output

true iff $x ∈ P$.

Algorithm

We check, using linear programming, the definition of a convex polytope that a point is in the set if and only if it is a convex combination of the vertices.

Let $\{v_j\}$ be the $m$ vertices of P. Then we solve the following $m$-dimensional linear program.

rand(::Type{VPolytope}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing,
     [num_vertices]::Int=-1)::VPolytope{N}

Create a random polytope in vertex representation.

Input

• VPolytope – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding
• num_vertices – (optional, default: -1) upper bound on the number of vertices of the polytope (see comment below)

Output

A random polytope in vertex representation.

Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1.

The number of vertices can be controlled with the argument num_vertices. For a negative value we choose a random number in the range dim:5*dim (except if dim == 1, in which case we choose in the range 1:2). Note that we do not guarantee that the vertices are not redundant.

translate(P::VPolytope{N}, v::AbstractVector{N}) where {N<:Real}

Translate (i.e., shift) a polytope in vertex representation by a given vector.

Input

• P – polytope in vertex representation
• v – translation vector

Output

A translated polytope in vertex representation.

Algorithm

We add the vector to each vertex of the polytope.

vertices_list(P::VPolytope{N})::Vector{Vector{N}} where {N<:Real}

Return the list of vertices of a polytope in V-representation.

Input

• P – polytope in vertex representation

Output

List of vertices.

remove_redundant_vertices(P::VPolytope{N};
                          [backend]=nothing,
                          [solver]=nothing)::VPolytope{N} where {N<:Real}

Return the polytope obtained by removing the redundant vertices of the given polytope.

Input

• P – polytope in vertex representation
• backend – (optional, default: nothing) the polyhedral computations backend, see default_polyhedra_backend(P1, N) or Polyhedra's documentation for further information on the available backends
• solver – (optional, default: nothing) the linear programming solver used in the backend, if needed; see default_lp_solver(N)

Output

A new polytope such that its vertices are the convex hull of the given polytope.

Notes

The optimization problem associated to removing redundant vertices is handled by Polyhedra. If the polyhedral computations backend requires an LP solver but it has not been set, we use default_lp_solver(N) to define such solver. Otherwise, the redundancy removal function of the polyhedral backend is used.

constraints_list(P::VPolytope{N}) where {N<:Real}

Return the list of constraints defining a polytope in V-representation.

Input

• P – polytope in V-representation

Output

The list of constraints of the polytope.

Algorithm

First the H-representation of $P$ is computed, then its list of constraints is returned.

tohrep(P::VPolytope{N};
       [backend]=default_polyhedra_backend(P, N)) where {N<:Real}

Transform a polytope in V-representation to a polytope in H-representation.

Input

• P – polytope in vertex representation
• backend – (optional, default: default_polyhedra_backend(P, N)) the backend for polyhedral computations

Output

The HPolytope which is the constraint representation of the given polytope in vertex representation.

Notes

The conversion may not preserve the numeric type (e.g., with N == Float32) depending on the backend. For further information on the supported backends see Polyhedra's documentation.

tovrep(P::VPolytope)

Return a vertex representation of the given polytope in vertex representation (no-op).

Input

• P – polytope in vertex representation

Output

The same polytope instance.

polyhedron(P::VPolytope{N};
           [backend]=default_polyhedra_backend(P, N)) where {N<:Real}

Return an VRep polyhedron from Polyhedra.jl given a polytope in V-representation.

Input

• P – polytope
• backend – (optional, default: default_polyhedra_backend(P, N)) the polyhedral computations backend, see Polyhedra's documentation for further information
• relative_dimension – (default, optional: nothing) an integer representing the (relative) dimension of the polytope; this argument is mandatory if the polytope is empty

Output

A VRep polyhedron.

Notes

The relative dimension (or just dimension) refers to the dimension of the set relative to itself, independently of the ambient dimension. For example, a point has (relative) dimension zero, and a line segment has (relative) dimension one.

In this library, LazySets.dim always returns the ambient dimension of the set, such that a line segment in two dimensions has dimension two. However, Polyhedra.dim will assign a dimension equal to one to a line segment because it uses a different convention.

linear_map(M::AbstractMatrix{N}, P::VPolytope{N}) where {N<:Real}

Concrete linear map of a polytope in vertex representation.

Input

• M – matrix
• P – polytope in vertex representation

Output

A polytope in vertex representation.

Algorithm

The linear map $M$ is applied to each vertex of the given set $P$, obtaining a polytope in V-representation. The output type is again a VPolytope.

PolynomialZonotope{N}

Type that represents a polynomial zonotope.

Fields

• c – starting point
• E – matrix of multi-indexed generators such that all indices have the same value
• F – matrix of multi-indexed generators such that not all indices have the same value
• G – matrix of single-indexed generators

Notes

Polynomial zonotopes were introduced by M. Althoff in [1] and have been applied as a non-convex set representation in the reachability problem of nonlinear ODEs.

Mathematically, a polynomial zonotope is defined as the tuple $(c, E, F, G)$, where:

• $c ∈ \mathbb{R}^n$ is the starting point (in some particular cases it corresponds to the center of a usual zonotope),

• $E = [E^{[1]} ⋯ E^{[η]}]$ is an $n × p × η(η+1)/2$ matrix with column-blocks

called the matrix of multi-indexed generators with equal indices, where each $f^{([i], k_1, k_2, …, k_i)}$ is an $n$-vector,

• $F = [F^{[2]} ⋯ F^{[η]}]$ is a matrix with column-blocks

called the matrix of multi-indexed generators with unequal indices (or, more accurately, not-all-equal indices), where each $f^{([i], k_1, k_2, …, k_i)}$ is an $n$-vector,

• $G = [G^{[1]} ⋯ G^{[q]}]$ is an $n × q$ matrix with columns

called the matrix of single-indexed generators, where each $g^{(i)}$ is an $n$-vector.

The polynomial zonotope $(c, E, F, G)$ defines the set:

where the number of factors in the final product, $β_j β_k ⋯ β_m$, corresponds to the polynomial order $η$.

[1] M. Althoff in Reachability analysis of nonlinear systems using conservative polynomialization and non-convex sets, Hybrid Systems: Computation and Control, 2013, pp. 173–182.

dim(pz::PolynomialZonotope)::Int

Return the ambient dimension of a polynomial zonotope.

Input

• pz – polynomial zonotope

Output

An integer representing the ambient dimension of the polynomial zonotope.

σ(d::AbstractVector{N}, pz::PolynomialZonotope{N})::AbstractVector{N} where {N}

Return the support vector of a polynomial zonotope along direction d.

Input

• d – direction
• pz – polynomial zonotope

Output

Vector representing the support vector.

ρ(d::AbstractVector{N}, pz::PolynomialZonotope{N})::AbstractVector{N} where {N}

Return the support function of a polynomial zonotope along direction d.

Input

• d – direction
• pz – polynomial zonotope

Output

Value of the support function.

polynomial_order(pz::PolynomialZonotope)::Int

Polynomial order of a polynomial zonotope.

Input

• pz – polynomial zonotope

Output

The polynomial order, defined as the maximal power of the scale factors $β_i$. Usually denoted $η$.

order(pz::PolynomialZonotope)::Rational{Int}

Order of a polynomial zonotope.

Input

• pz – polynomial zonotope

Output

The order, a rational number defined as the total number of generators divided by the ambient dimension.

linear_map(M::Matrix, pz::PolynomialZonotope)

Return the linear map of a polynomial zonotope.

Input

• M – matrix
• pz – polynomial zonotope

Output

Polynomial zonotope such that its starting point and generators are those of pz multiplied by the matrix M.

scale(α::Number, pz::PolynomialZonotope)

Return a polynomial zonotope modified by a scale factor.

Input

• α – polynomial zonotope
• pz – polynomial zonotope

Output

Polynomial zonotope such that its center and generators are multiples of those of pz by a factor $α$.

Singleton{N<:Real, VN<:AbstractVector{N}} <: AbstractSingleton{N}

Type that represents a singleton, that is, a set with a unique element.

Fields

• element – the only element of the set

rand(::Type{Singleton}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::Singleton{N}

Create a random singleton.

Input

• Singleton – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A random singleton.

Algorithm

The element is a normally distributed vector with entries of mean 0 and standard deviation 1.

element(S::Singleton{N}) where {N<:Real}

Return the element of a singleton.

Input

• S – singleton

Output

The element of the singleton.

element(S::Singleton{N}, i::Int)::N where {N<:Real}

Return the i-th entry of the element of a singleton.

Input

• S – singleton
• i – dimension

Output

The i-th entry of the element of the singleton.

translate(S::Singleton{N}, v::AbstractVector{N}) where {N<:Real}

Translate (i.e., shift) a singleton by a given vector.

Input

• S – singleton
• v – translation vector

Output

A translated singleton.

Algorithm

We add the vector to the point in the singleton.

Universe{N<:Real} <: AbstractPolyhedron{N}

Type that represents the universal set, i.e., the set of all elements.

dim(U::Universe)

Return the dimension of a universe.

Input

• U – universe

Output

The dimension of a universe.

ρ(d::AbstractVector{N}, U::Universe{N}) where {N<:Real}

Return the support function of a universe.

Input

• d – direction
• U – universe

Output

The support function in the given direction.

Algorithm

If the direction is all zero, the result is zero. Otherwise, the result is Inf.

σ(d::AbstractVector{N}, U::Universe{N}) where {N<:Real}

Return the support vector of a universe.

Input

• d – direction
• U – universe

Output

A vector with infinity values, except in dimensions where the direction is zero.

∈(x::AbstractVector{N}, U::Universe{N})::Bool where {N<:Real}

Check whether a given point is contained in a universe.

Input

• x – point/vector
• U – universe

Output

The output is always true.

Examples

julia> [1.0, 0.0] ∈ Universe(2)
true

rand(::Type{Universe}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::Universe{N}

Create a universe (note that there is nothing to randomize).

Input

• Universe – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

The (only) universe of the given numeric type and dimension.

an_element(U::Universe{N}) where {N<:Real}

Return some element of a universe.

Input

• U – universe

Output

The origin.

isempty(U::Universe)::Bool

Return if a universe is empty or not.

Input

• U – universe

Output

false.

isbounded(U::Universe)::Bool

Determine whether a universe is bounded.

Input

• U – universe

Output

false as the universe is unbounded.

isuniversal(U::Universe{N}, [witness]::Bool=false
           )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether a universe is universal.

Input

• U – universe
• witness – (optional, default: false) compute a witness if activated

Output

• If witness option is deactivated: true
• If witness option is activated: (true, [])

norm(U::Universe, [p]::Real=Inf)

Return the norm of a universe. It is the norm of the enclosing ball (of the given $p$-norm) of minimal volume that is centered in the origin.

Input

• U – universe
• p – (optional, default: Inf) norm

Output

An error.

radius(U::Universe, [p]::Real=Inf)

Return the radius of a universe. It is the radius of the enclosing ball (of the given $p$-norm) of minimal volume with the same center.

Input

• U – universe
• p – (optional, default: Inf) norm

Output

An error.

diameter(U::Universe, [p]::Real=Inf)

Return the diameter of a universe. It is the maximum distance between any two elements of the set, or, equivalently, the diameter of the enclosing ball (of the given $p$-norm) of minimal volume with the same center.

Input

• U – universe
• p – (optional, default: Inf) norm

Output

An error.

constraints_list(U::Universe{N}) where {N<:Real}

Return the list of constraints defining a universe.

Input

• U – universe

Output

The empty list of constraints, as the universe is unconstrained.

constrained_dimensions(U::Universe)::Vector{Int}

Return the indices in which a universe is constrained.

Input

• U – universe

Output

The empty vector, as the universe is unconstrained in every dimension.

translate(U::Universe{N}, v::AbstractVector{N}) where {N<:Real}

Translate (i.e., shift) a universe by a given vector.

Input

• U – universe
• v – translation vector

Output

The universe.

ZeroSet{N<:Real} <: AbstractSingleton{N}

Type that represents the zero set, i.e., the set that only contains the origin.

Fields

• dim – the ambient dimension of this zero set

dim(Z::ZeroSet)::Int

Return the ambient dimension of this zero set.

Input

• Z – a zero set, i.e., a set that only contains the origin

Output

The ambient dimension of the zero set.

σ(d::AbstractVector{N}, Z::ZeroSet{N}) where {N<:Real}

Return the support vector of a zero set.

Input

• d – direction
• Z – a zero set, i.e., a set that only contains the origin

Output

The returned value is the origin since it is the only point that belongs to this set.

ρ(d::AbstractVector{N}, Z::ZeroSet{N}) where {N<:Real}

Evaluate the support function of a zero set in a given direction.

Input

• d – direction
• Z – a zero set, i.e., a set that only contains the origin

Output

0.

∈(x::AbstractVector{N}, Z::ZeroSet{N})::Bool where {N<:Real}

Check whether a given point is contained in a zero set.

Input

• x – point/vector
• Z – zero set

Output

true iff $x ∈ Z$.

Examples

julia> Z = ZeroSet(2);

julia> [1.0, 0.0] ∈ Z
false
julia> [0.0, 0.0] ∈ Z
true

rand(::Type{ZeroSet}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::ZeroSet{N}

Create a zero set (note that there is nothing to randomize).

Input

• ZeroSet – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

The (only) zero set of the given numeric type and dimension.

element(S::ZeroSet{N})::Vector{N} where {N<:Real}

Return the element of a zero set.

Input

• S – zero set

Output

The element of the zero set, i.e., a zero vector.

element(S::ZeroSet{N}, ::Int)::N where {N<:Real}

Return the i-th entry of the element of a zero set.

Input

• S – zero set
• i – dimension

Output

The i-th entry of the element of the zero set, i.e., 0.

linear_map(M::AbstractMatrix{N}, Z::ZeroSet{N}) where {N<:Real}

Concrete linear map of a zero set.

Input

• M – matrix
• Z – zero set

Output

The zero set whose dimension matches the output dimension of the given matrix.

translate(Z::ZeroSet{N}, v::AbstractVector{N}) where {N<:Real}

Translate (i.e., shift) a zero set by a given vector.

Input

• Z – zero set
• v – translation vector

Output

A singleton containing the vector v.

Zonotope{N<:Real} <: AbstractZonotope{N}

Type that represents a zonotope.

Fields

• center – center of the zonotope
• generators – matrix; each column is a generator of the zonotope

Notes

Mathematically, a zonotope is defined as the set

where $c \in \mathbb{R}^n$ is its center and $\{g_i\}_{i=1}^p$, $g_i \in \mathbb{R}^n$, is the set of generators. This characterization defines a zonotope as the finite Minkowski sum of line segments. Zonotopes can be equivalently described as the image of a unit infinity-norm ball in $\mathbb{R}^n$ by an affine transformation.

Zonotopes can be constructed in two different ways: either passing the generators as a matrix, where each column represents a generator, or passing a list of vectors where each vector represents a generator. Below we illustrate both ways.

The optional argument remove_zero_generators controls whether we remove zero columns from the generators matrix. This option is active by default.

Examples

A two-dimensional zonotope with given center and set of generators:

julia> Z = Zonotope([1.0, 0.0], [0.1 0.0; 0.0 0.1])
Zonotope{Float64}([1.0, 0.0], [0.1 0.0; 0.0 0.1])

julia> dim(Z)
2

Here, each column of the second input corresponds to a generator.

We can collect its vertices using vertices_list:

julia> vertices_list(Z)
4-element Array{Array{Float64,1},1}:
 [1.1, 0.1]
 [0.9, 0.1]
 [0.9, -0.1]
 [1.1, -0.1] 

The support vector along a given direction can be computed using σ (resp. the support function can be computed using ρ):

julia> σ([1., 1.], Z)
2-element Array{Float64,1}:
 1.1
 0.1

Zonotopes admit an alternative constructor that receives a list of vectors, each vector representing a generator:

julia> Z = Zonotope(ones(2), [[1., 0.], [0., 1.], [1., 1.]])
Zonotope{Float64}([1.0, 1.0], [1.0 0.0 1.0; 0.0 1.0 1.0])

julia> Z.generators
2×3 Array{Float64,2}:
 1.0  0.0  1.0
 0.0  1.0  1.0

ρ(d::AbstractVector{N}, Z::Zonotope{N}) where {N<:Real}

Return the support function of a zonotope in a given direction.

Input

• d – direction
• Z – zonotope

Output

The support function of the zonotope in the given direction.

Algorithm

The support value is $cᵀ d + ‖Gᵀ d‖₁$ where $c$ is the center and $G$ is the generator matrix of Z.

σ(d::AbstractVector{N}, Z::Zonotope{N}) where {N<:Real}

Return the support vector of a zonotope in a given direction.

Input

• d – direction
• Z – zonotope

Output

Support vector in the given direction. If the direction has norm zero, the vertex with $ξ_i = 1 \ \ ∀ i = 1,…, p$ is returned.

∈(x::AbstractVector{N}, Z::Zonotope{N};
  solver=GLPKSolverLP(method=:Simplex))::Bool where {N<:Real}

Check whether a given point is contained in a zonotope.

Input

• x – point/vector
• Z – zonotope
• solver – (optional, default: GLPKSolverLP(method=:Simplex)) the backend used to solve the linear program

Output

true iff $x ∈ Z$.

Examples

julia> Z = Zonotope([1.0, 0.0], [0.1 0.0; 0.0 0.1]);

julia> [1.0, 0.2] ∈ Z
false
julia> [1.0, 0.1] ∈ Z
true

Algorithm

The membership problem is computed by stating and solving the following linear program with the simplex method. Let $p$ and $n$ be the number of generators and ambient dimension, respectively. We consider the minimization of $x_0$ in the $p+1$-dimensional space of elements $(x_0, ξ_1, …, ξ_p)$ constrained to $0 ≤ x_0 ≤ ∞$, $ξ_i ∈ [-1, 1]$ for all $i = 1, …, p$, and such that $x-c = Gξ$ holds. If a feasible solution exists, the optimal value $x_0 = 0$ is achieved.

Notes

This function is parametric in the number type N. For exact arithmetic use an appropriate backend, e.g. solver=GLPKSolverLP(method=:Exact).

rand(::Type{Zonotope}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
     [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
    )::Zonotope{N}

Create a random zonotope.

Input

• Zonotope – type for dispatch
• N – (optional, default: Float64) numeric type
• dim – (optional, default: 2) dimension
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding
• num_generators – (optional, default: -1) number of generators of the zonotope (see comment below)

Output

A random zonotope.

Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1.

The number of generators can be controlled with the argument num_generators. For a negative value we choose a random number in the range dim:2*dim (except if dim == 1, in which case we only create a single generator).

vertices_list(Z::Zonotope{N};
              [apply_convex_hull]::Bool=true)::Vector{Vector{N}} where {N<:Real}

Return the vertices of a zonotope.

Input

• Z – zonotope
• apply_convex_hull – (optional, default: true) if true, post-process the computation with the convex hull of the points

Output

List of vertices as a vector of vectors.

Algorithm

If the zonotope has $p$ generators, each vertex is the result of summing the center with some linear combination of generators, where the combination factors are $ξ_i ∈ \{-1, 1\}$.

There are at most $2^p$ distinct vertices. Use the flag apply_convex_hull to control whether a convex hull algorithm is applied to the vertices computed by this method; otherwise, redundant vertices may be present.

constraints_list(P::Zonotope{N}) where {N<:Real}

Return the list of constraints defining a zonotope.

Input

• Z – zonotope

Output

The list of constraints of the zonotope.

Algorithm

This is the (inefficient) fallback implementation for rational numbers. It first computes the vertices and then converts the corresponding polytope to constraint representation.

constraints_list(Z::Zonotope{N}) where {N<:AbstractFloat}

Return the list of constraints defining a zonotope.

Input

• Z – zonotope
• check_full_rank – (optional; default: true) flag for checking whether the generator matrix has full rank

Output

The list of constraints of the zonotope.

Notes

The algorithm assumes that no generator is redundant. The result has $2 \binom{p}{n-1}$ (with $p$ being the number of generators and $n$ being the ambient dimension) constraints, which is optimal under this assumption.

If $p < n$ or the generator matrix is not full rank, we fall back to the (slower) computation based on the vertex representation.

Algorithm

We follow the algorithm presented in Althoff, Stursberg, Buss: Computing Reachable Sets of Hybrid Systems Using a Combination of Zonotopes and Polytopes. 2009.

The one-dimensional case is not covered by that algorithm; we manually handle this case, assuming that there is only one generator.

center(Z::Zonotope{N})::Vector{N} where {N<:Real}

Return the center of a zonotope.

Input

• Z – zonotope

Output

The center of the zonotope.

order(Z::Zonotope)::Rational

Return the order of a zonotope.

Input

• Z – zonotope

Output

A rational number representing the order of the zonotope.

Notes

The order of a zonotope is defined as the quotient of its number of generators and its dimension.

generators(Z::Zonotope)

Return an iterator over the generators of a zonotope.

Input

• Z – zonotope

Output

An iterator over the generators of Z.

genmat(Z::Zonotope)

Return the generator matrix of a zonotope.

Input

• Z – zonotope

Output

A matrix where each column represents one generator of the zonotope Z.

linear_map(M::AbstractMatrix{N}, Z::Zonotope{N}) where {N<:Real}

Concrete linear map of a zonotope.

Input

• M – matrix
• Z – zonotope

Output

The zonotope obtained by applying the linear map to the center and generators of $Z$.

translate(Z::Zonotope{N}, v::AbstractVector{N}; share::Bool=false
         ) where {N<:Real}

Translate (i.e., shift) a zonotope by a given vector.

Input

• Z – zonotope
• v – translation vector
• share – (optional, default: false) flag for sharing unmodified parts of the original set representation

Output

A translated zonotope.

Notes

The generator matrix is shared with the original zonotope if share == true.

Algorithm

We add the vector to the center of the zonotope.

scale(α::Real, Z::Zonotope)

Concrete scaling of a zonotope.

Input

• α – scalar
• Z – zonotope

Output

The zonotope obtained by applying the numerical scale to the center and generators of $Z$.

ngens(Z::Zonotope)::Int

Return the number of generators of a zonotope.

Input

• Z – zonotope

Output

Integer representing the number of generators.

reduce_order(Z::Zonotope, r)::Zonotope

Reduce the order of a zonotope by overapproximating with a zonotope with less generators.

Input

• Z – zonotope
• r – desired order

Output

A new zonotope with less generators, if possible.

Algorithm

This function implements the algorithm described in A. Girard's Reachability of Uncertain Linear Systems Using Zonotopes, HSCC. Vol. 5. 2005.

If the desired order is smaller than one, the zonotope is not reduced.

split(Z::Zonotope, j::Int)

Return two zonotopes obtained by splitting the given zonotope.

Input

• Z – zonotope
• j – index of the generator to be split

Output

The zonotope obtained by splitting Z into two zonotopes such that their union is Z and their intersection is possibly non-empty.

Algorithm

This function implements [Prop. 3, 1], that we state next. The zonotope $Z = ⟨c, g^{(1, …, p)}⟩$ is split into:

such that $Z₁ ∪ Z₂ = Z$ and $Z₁ ∩ Z₂ = Z^*$, where

[1] Althoff, M., Stursberg, O., & Buss, M. (2008). Reachability analysis of nonlinear systems with uncertain parameters using conservative linearization. In Proc. of the 47th IEEE Conference on Decision and Control.