Approximations

This section of the manual describes the Cartesian decomposition algorithms and the approximation of high-dimensional convex sets using projections.

Cartesian Decomposition

LazySets.Approximations.decomposeFunction
decompose(S::ConvexSet{N},
          partition::AbstractVector{<:AbstractVector{Int}},
          block_options
         ) where {N}

Decompose a high-dimensional set into a Cartesian product of overapproximations of the projections over the specified subspaces.

Input

  • S – set
  • partition – vector of blocks (i.e., of vectors of integers) (see the Notes below)
  • block_options – mapping from block indices in partition to a corresponding overapproximation option; we only require access via [⋅] (but see also the Notes below)

Output

A CartesianProductArray containing the low-dimensional approximated projections.

Algorithm

For each block a specific project method is called, dispatching on the corresponding overapproximation option.

Notes

The argument partition requires some discussion. Typically, the list of blocks should form a partition of the set $\{1, \dots, n\}$ represented as a list of consecutive blocks, where $n$ is the ambient dimension of set S.

However, technically there is no problem if the blocks are not consecutive, blocks are missing, blocks occur more than once, or blocks are overlapping. This function will, however, stick to the order of blocks, so the resulting set must be interpreted with care in such cases. One use case is the need of a projection consisting of several blocks.

For convenience, the argument block_options can also be given as a single option instead of a mapping, which is then interpreted as the option for all blocks.

Examples

This function supports different options: one can specify the target set, the degree of accuracy, and template directions. These options are exemplified below, where we use the following example.

julia> S = Ball2(zeros(4), 1.0);  # set to be decomposed (4D 2-norm unit ball)

julia> P2d = [1:2, 3:4];  # a partition with two blocks of size two

julia> P1d = [[1], [2], [3], [4]];  # a partition with four blocks of size one

Different set types

We can decompose using polygons in constraint representation:

julia> all(ai isa HPolygon for ai in array(decompose(S, P2d, HPolygon)))
true

For decomposition into 1D subspaces, we can use Interval:

julia> all(ai isa Interval for ai in array(decompose(S, P1d, Interval)))
true

However, if you need to specify different set types for different blocks, the interface presented so far does not apply. See the paragraph Advanced input for different block approximations below for how to do that.

Refining the decomposition I: $ε$-close approximation

The $ε$ option can be used to refine a decomposition, i.e., obtain a more accurate result. We use the Iterative refinement algorithm from the Approximations module.

To illustrate this, consider again the set S from above. We decompose into two 2D polygons. Using smaller $ε$ implies a better precision, thus more constraints in each 2D decomposition. In the following example, we look at the number of constraints in the first block.

julia> d(ε, bi) = array(decompose(S, P2d, (HPolygon => ε)))[bi]
d (generic function with 1 method)

julia> [length(constraints_list(d(ε, 1))) for ε in [Inf, 0.1, 0.01]]
3-element Vector{Int64}:
  4
  8
 32

Refining the decomposition II: template polyhedra

Another way to refine a decomposition is by using template polyhedra. The idea is to specify a set of template directions and then to compute on each block the polytopic overapproximation obtained by evaluating the support function of the given input set over the template directions.

For example, octagonal 2D approximations of the set S are obtained with:

julia> B = decompose(S, P2d, OctDirections);

julia> length(B.array) == 2 && all(dim(bi) == 2 for bi in B.array)
true

See Template directions for the available template directions. Note that, in contrast to the polygonal $ε$-close approximation from above, this method can be applied to blocks of any size.

julia> B = decompose(S, [1:4], OctDirections);

julia> length(B.array) == 1 && dim(B.array[1]) == 4
true

Advanced input for different block approximations

Instead of defining the approximation option uniformly for each block, we can define different approximations for different blocks. The third argument has to be a mapping from block index (in the partition) to the corresponding approximation option.

For example:

julia> res = array(decompose(S, P2d, Dict(1 => Hyperrectangle, 2 => 0.1)));

julia> typeof(res[1]), typeof(res[2])
(Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}, HPolygon{Float64, Vector{Float64}})
source

Convenience functions

LazySets.Approximations.uniform_partitionFunction
 uniform_partition(n::Int, block_size::Int)

Compute a uniform block partition of the given size.

Input

  • n – number of dimensions of the partition
  • block_size – size of each block

Output

A vector of ranges, Vector{UnitRange{Int}}, such that the size of each block is the same, if possible.

Examples

If the number of dimensions n is 2, we have two options: either two blocks of size 1 or one block of size 2:

julia> LazySets.Approximations.uniform_partition(2, 1)
2-element Vector{UnitRange{Int64}}:
 1:1
 2:2

julia> LazySets.Approximations.uniform_partition(2, 2)
1-element Vector{UnitRange{Int64}}:
 1:2

If the block size argument is not compatible with (i.e. does not divide) n, the output is filled with one block of the size needed to reach n:

julia> LazySets.Approximations.uniform_partition(3, 1)
3-element Vector{UnitRange{Int64}}:
 1:1
 2:2
 3:3

julia> LazySets.Approximations.uniform_partition(3, 2)
2-element Vector{UnitRange{Int64}}:
 1:2
 3:3

julia> LazySets.Approximations.uniform_partition(10, 6)
2-element Vector{UnitRange{Int64}}:
 1:6
 7:10
source

Overapproximations

LazySets.Approximations.overapproximateFunction
overapproximate(X::S, ::Type{S}, args...) where {S<:ConvexSet}

Overapproximating a set of type S with type S is a no-op.

Input

  • X – set
  • Type{S} – target set type
  • args – further arguments (ignored)

Output

The input set.

source
overapproximate(S::ConvexSet)

Alias for overapproximate(S, Hyperrectangle) resp. box_approximation(S).

source
overapproximate(S::ConvexSet, ::Type{<:Hyperrectangle})

Alias for box_approximation(S).

source
overapproximate(S::ConvexSet, ::Type{<:BallInf})

Alias for ballinf_approximation(S).

source
overapproximate(S::ConvexSet{N},
                ::Type{<:HPolygon},
                [ε]::Real=Inf) where {N}

Return an approximation of a given 2D set using iterative refinement.

Input

  • S – convex set, assumed to be two-dimensional
  • HPolygon – type for dispatch
  • ε – (optional, default: Inf) error tolerance
  • prune – (optional, default: true) flag for removing redundant constraints in the end

Output

A polygon in constraint representation.

Notes

The result is always a convex overapproximation of the input set.

If no error tolerance ε is given, or is Inf, the result is a box-shaped polygon. For convex input sets, the result is an ε-close approximation as a polygon, with respect to the Hausdorff distance.

source
overapproximate(S::ConvexSet, ε::Real)

Alias for overapproximate(S, HPolygon, ε).

source
overapproximate(X::ConvexHull{N, <:AbstractZonotope, <:AbstractZonotope},
                ::Type{<:Zonotope}) where {N}

Overapproximate the convex hull of two zonotopes.

Input

  • X – convex hull of two zonotopes
  • Zonotope – type for dispatch
  • algorithm – (optional; default: "mean") choice of algorithm; possible values are "mean" and "join"

Output

A zonotope $Z$ such that $X ⊆ Z$.

Algorithm

The algorithm can be controlled by the parameter algorithm. Note that the results of the two implemented algorithms are generally incomparable.

'mean' method

If algorithm == "mean", we choose the method proposed in [1]. The convex hull of two zonotopes $Z₁$ and $Z₂$ of the same order, which we write

\[Z_j = ⟨c^{(j)}, g^{(j)}_1, …, g^{(j)}_p⟩\]

for $j = 1, 2$, can be overapproximated as follows:

\[CH(Z_1, Z_2) ⊆ \frac{1}{2}⟨c^{(1)}+c^{(2)}, g^{(1)}_1+g^{(2)}_1, …, g^{(1)}_p+g^{(2)}_p, c^{(1)}-c^{(2)}, g^{(1)}_1-g^{(2)}_1, …, g^{(1)}_p-g^{(2)}_p⟩.\]

If the zonotope order is not the same, this algorithm calls reduce_order to reduce the order to the minimum of the arguments.

It should be noted that the output zonotope is not necessarily the minimal enclosing zonotope, which is in general expensive in high dimensions. This is further investigated in [2].

'join' method

If algorithm == "join", we choose the method proposed in [3, Definition 1]. The convex hull $X$ of two zonotopes $Z₁$ and $Z₂$ is overapproximated by a zonotope $Z₃$ such that the box approximation of $X$ is identical with the box approximation of $Z₃$. Let $□(X)$ denote the box approximation of $X$. The center of $Z₃$ is the center of $□(X)$.

The generator construction consists of two phases. In the first phase, we construct generators $g$ as a combination of one generator from $Z₁$, say, $g₁$, with another generator from $Z₂$, say, $g₂$. The entry of $g$ in the $i$-th dimension is given as

\[ g[i] = \arg\min_{\min(g₁[i], g₂[i]) ≤ x ≤ \max(g₁[i], g₂[i])} |x|.\]

If $g$ is the zero vector, it can be omitted.

In the second phase, we construct another generator for each dimension. These generators are scaled unit vectors. The following formula defines the sum of all those generators.

\[ \sup(□(X)) - c - ∑_g |g|\]

where $c$ is the center of the new zonotope and the $g$s are the generators constructed in the first phase.

References

[1] Reachability of Uncertain Linear Systems Using Zonotopes, A. Girard. HSCC 2005.

[2] Zonotopes as bounding volumes, L. J. Guibas et al, Proc. of Symposium on Discrete Algorithms, pp. 803-812.

[3] The zonotope abstract domain Taylor1+. K. Ghorbal, E. Goubault, S. Putot. CAV 2009.

source
overapproximate(lm::LinearMap{N, <:AbstractZonotope},
                ::Type{<:Zonotope}) where {N}

Overapproximate a lazy linear map of a zonotopic set with a zonotope.

Input

  • lm – lazy linear map of a zonotopic set
  • Zonotope – type for dispatch

Output

The tight zonotope corresponding to lm.

source
overapproximate(X::ConvexSet{N}, dir::AbstractDirections; [prune]::Bool=true) where {N}

Overapproximate a (possibly unbounded) set with template directions.

Input

  • X – set
  • dir – (concrete) direction representation
  • prune – (optional, default: true) flag for removing redundant constraints in the end

Output

A polyhedron overapproximating the set X with the directions from dir. The overapproximation is computed using support functions. If the obtained set is bounded, the result is an HPolytope. Otherwise the result is an HPolyhedron.

source
overapproximate(X::ConvexSet{N}, dir::Type{<:AbstractDirections}) where {N}

Overapproximating a set with template directions.

Input

  • X – set
  • dir – type of direction representation

Output

A polyhedron overapproximating the set X with the directions from dir. If the directions are known to be bounded, the result is an HPolytope, otherwise the result is an HPolyhedron.

source
overapproximate(S::ConvexSet{N}, ::Type{<:Interval}) where {N}

Return the overapproximation of a unidimensional set with an interval.

Input

  • S – one-dimensional set
  • Interval – type for dispatch

Output

An interval.

Algorithm

We use two support-function evaluations.

source
overapproximate(cap::Intersection, ::Type{<:Interval})

Return the overapproximation of a unidimensional intersection with an interval.

Input

  • cap – one-dimensional lazy intersection
  • Interval – type for dispatch

Output

An interval.

Algorithm

The algorithm recursively overapproximates the two intersected sets with intervals and then intersects these.

source
overapproximate(cap::IntersectionArray, ::Type{<:Interval})

Return the overapproximation of a unidimensional intersection with an interval.

Input

  • cap – one-dimensional lazy intersection
  • Interval – type for dispatch

Output

An interval.

Algorithm

The algorithm recursively overapproximates the two intersected sets with intervals and then intersects these.

source
overapproximate(cap::Intersection{N, <:ConvexSet, <:AbstractPolyhedron},
                dir::AbstractDirections;
                kwargs...
               ) where {N}

Return the overapproximation of the intersection between a compact set and a polytope given a set of template directions.

Input

  • cap – intersection of a compact set and a polytope
  • dir – template directions
  • kwargs – additional arguments that are passed to the support function algorithm

Output

A polytope in H-representation such that the normal direction of each half-space is given by an element of dir.

Algorithm

Let di be a direction drawn from the set of template directions dir. Let X be the compact set and let P be the polytope. We overapproximate the set X ∩ H with a polytope in constraint representation using a given set of template directions dir.

The idea is to solve the univariate optimization problem ρ(di, X ∩ Hi) for each half-space in the set P and then take the minimum. This gives an overapproximation of the exact support function.

This algorithm is inspired from G. Frehse, R. Ray. Flowpipe-Guard Intersection for Reachability Computations with Support Functions.

Notes

This method relies on having available the constraints_list of the polytope P.

This method of overapproximations can return a non-empty set even if the original intersection is empty.

source
overapproximate(cap::Intersection{N, <:HalfSpace, <:AbstractPolytope},
                dir::AbstractDirections;
                [kwargs]...
               ) where {N}

Return the overapproximation of the intersection between a half-space and a polytope given a set of template directions.

Input

  • cap – intersection of a half-space and a polytope
  • dir – template directions
  • kwargs – additional arguments that are passed to the support function algorithm

Output

A polytope in H-representation such that the normal direction of each half-space is given by an element of dir.

source
overapproximate(P::SimpleSparsePolynomialZonotope, ::Type{Zonotope}; nsdiv=1, partition=nothing)

Return a zonotope containing $P$.

Input

  • P – simple sparse polynomial zonotope
  • nsdiv – (optional, default: 1) size of uniform partitioning grid
  • partition – (optional, default: nothing) tuple of integers indicating the number of partitions in each dimensino, the length should match nparams(P)

Output

A zonotope containing P.

source
overapproximate(P::SimpleSparsePolynomialZonotope, ::Type{Zonotope}, dom::IntervalBox)

Compute the zonotope overapproximation of the given sparse polynomial zonotope over the parameter domain dom, which should be a subset of [-1, 1]^q, where q = nparams(P).

source
overapproximate(lm::LinearMap{N, <:CartesianProductArray},
                ::Type{CartesianProductArray{N, S}}
               ) where {N, S<:ConvexSet}

Decompose a lazy linear map of a Cartesian product array while keeping the original block structure.

Input

  • lm – lazy linear map of Cartesian product array
  • CartesianProductArray – type for dispatch

Output

A CartesianProductArray representing the decomposed linear map.

source
overapproximate(lm::LinearMap{N, <:CartesianProductArray},
                ::Type{<:CartesianProductArray},
                dir::Type{<:AbstractDirections}) where {N}

Decompose a lazy linear map of a Cartesian product array with template directions while keeping the original block structure.

Input

  • lm – lazy linear map of a Cartesian product array
  • CartesianProductArray – type for dispatch
  • dir – template directions for overapproximation

Output

A CartesianProductArray representing the decomposed linear map.

source
overapproximate(lm::LinearMap{N, <:CartesianProductArray},
                ::Type{<:CartesianProductArray},
                set_type::Type{<:ConvexSet}) where {N}

Decompose a lazy linear map of a Cartesian product array with a given set type while keeping the original block structure.

Input

  • lm – lazy linear map of a Cartesian product array
  • CartesianProductArray – type for dispatch
  • set_type – set type for overapproximation

Output

A CartesianProductArray representing the decomposed linear map.

source
overapproximate(rm::ResetMap{N, <:CartesianProductArray},
                ::Type{<:CartesianProductArray}, oa) where {N}

Overapproximate a reset map (that only resets to zero) of a Cartesian product by a new Cartesian product.

Input

  • rm – reset map
  • CartesianProductArray – type for dispatch
  • oa – overapproximation option

Output

A Cartesian product with the same block structure.

Notes

This implementation currently only supports resets to zero.

Algorithm

We convert the ResetMap into a LinearMap and then call the corresponding overapproximation method.

source
overapproximate(cap::Intersection{N,
                                  <:CartesianProductArray,
                                  <:AbstractPolyhedron},
                ::Type{CartesianProductArray}, oa) where {N}

Return the intersection of the Cartesian product of a finite number of convex sets and a polyhedron.

Input

  • cap – lazy intersection of a Cartesian product array and a polyhedron
  • CartesianProductArray – type for dispatch
  • oa – overapproximation option

Output

A CartesianProductArray that overapproximates the intersection of cpa and P.

Algorithm

The intersection only needs to be computed in the blocks of cpa that are constrained in P. Hence we first collect those constrained blocks in a lower-dimensional Cartesian product array and then convert to an HPolytope X. Then we take the intersection of X and the projection of Y onto the corresponding dimensions. (This projection is purely syntactic and exact.) Finally we decompose the result again and plug together the unaffected old blocks and the newly computed blocks. The result is a CartesianProductArray with the same block structure as in X.

source
overapproximate(Z::Zonotope{N}, ::Type{<:Zonotope}, r::Union{Integer, Rational}) where {N}

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

Input

  • Z – zonotope
  • Zonotope – desired type for dispatch
  • r – desired order

Output

A new zonotope with less generators, if possible.

Algorithm

This function falls back to reduce_order with the default algorithm.

source
overapproximate(X::ConvexSet, ZT::Type{<:Zonotope},
                dir::AbstractDirections;
                [algorithm]="vrep", kwargs...)

Overapproximate a polytopic set with a zonotope.

Input

  • X – polytopic set
  • Zonotope – type for dispatch
  • dir – directions used for the generators
  • algorithm – (optional, default: "vrep") method used to compute the overapproximation
  • kwargs – further algorithm choices

Output

A zonotope that overapproximates X and uses at most the directions provided in dir (redundant directions will be ignored).

Notes

Two algorithms are available:

  • "vrep" – Overapproximate a polytopic set with a zonotope of minimal total generator sum using only generators in the given directions. Under this constraint, the zonotope has the minimal sum of generator vectors. See the docstring of _overapproximate_zonotope_vrep for further details.

  • "cpa" – Overapproximate a polytopic set with a zonotope using a cartesian decomposition into two-dimensional blocks. See the docstring of _overapproximate_zonotope_cpa for further details.

source
overapproximate(r::Rectification{N, <:AbstractZonotope}, ::Type{<:Zonotope}) where {N}

Overapproximation of the rectification of a zonotopic set.

Input

  • r – lazy rectification of a zonotopic set
  • Zonotope – type for dispatch

Output

A zonotope overapproximation of the set obtained by rectifying Z.

Algorithm

This function implements [Theorem 3.1, 1].

[1] Singh, G., Gehr, T., Mirman, M., Püschel, M., & Vechev, M. (2018). Fast and effective robustness certification. In Advances in Neural Information Processing Systems (pp. 10802-10813).

source
overapproximate(CHA::ConvexHullArray{N, <:AbstractZonotope}, ::Type{<:Zonotope}) where {N}

Overapproximation of the convex hull array of zonotopic sets.

Input

  • CHA – convex hull array of zonotopic sets
  • Zonotope – type for dispatch

Output

A zonotope overapproximation of the convex hull array of zonotopic sets.

Algorithm

This function iteratively applies the overapproximation algorithm for the convex hull of two zonotopes to the given array of zonotopes.

source
overapproximate(Z::AbstractZonotope, ::Type{<:HParallelotope}, indices=1:dim(Z))

Overapproximation of a zonotopic set with a parallelotopic set in constraint representation.

Input

  • Z – zonotopic set
  • HParallelotope – type for dispatch
  • indices – (optional; default: 1:dim(Z)) generator indices selected when constructing the parallelotope

Output

An overapproximation of the given zonotope using a parallelotope.

Algorithm

The algorithm is based on Proposition 8 discussed in Section 5 of [1].

[1] Althoff, M., Stursberg, O., & Buss, M. (2010). Computing reachable sets of hybrid systems using a combination of zonotopes and polytopes. Nonlinear analysis: hybrid systems, 4(2), 233-249.

source
overapproximate(X::Intersection{N, <:AbstractZonotope, <:Hyperplane},
                dirs::AbstractDirections) where {N}

Overapproximation of the intersection between a zonotopic set and a hyperplane

Input

  • X – intersection between a zonotopic set and a hyperplane
  • dirs – type of direction representation

Output

An overapproximation of the intersection between a zonotopic set and a hyperplane.

Algorithm

This function implements [Algorithm 8.1, 1].

[1] Colas Le Guernic. Reachability Analysis of Hybrid Systems with Linear Continuous Dynamics. Computer Science [cs]. Université Joseph-Fourier - Grenoble I, 2009. English. fftel-00422569v2f

source
overapproximate(qm::QuadraticMap{N, <:AbstractZonotope}, ::Type{<:Zonotope}) where {N}

Return a zonotopic overapproximation of the quadratic map of the given zonotope.

Input

  • Z – zonotope
  • Q – array of square matrices

Output

A zonotopic overapproximation of the quadratic map of the given zonotope.

Notes

Mathematically, a quadratic map of a zonotope is defined as:

\[Z_Q = \right\{ \lambda | \lambda_i = x^T Q\^{(i)} x,~i = 1, \ldots, n,~x \in Z \left\}\]

Algorithm

This function implements [Lemma 1, 1].

[1] Matthias Althoff and Bruce H. Krogh. 2012. Avoiding geometric intersection operations in reachability analysis of hybrid systems. In Proceedings of the 15th ACM international conference on Hybrid Systems: Computation and Control (HSCC ’12). Association for Computing Machinery, New York, NY, USA, 45–54.

source
overapproximate(vTM::Vector{TaylorModel1{T, S}}, ::Type{<:Zonotope};
                [remove_zero_generators]::Bool=true
                [normalize]::Bool=true) where {T, S}

Overapproximate a taylor model in one variable with a zonotope.

Input

  • vTMTaylorModel1
  • Zonotope – type for dispatch
  • remove_zero_generators – (optional; default: true) flag to remove zero generators of the resulting zonotope
  • normalize – (optional; default: true) flag to skip the normalization of the Taylor models

Output

A zonotope that overapproximates the range of the given taylor model.

Examples

If the polynomials are linear, this functions exactly transforms to a zonotope. However, the nonlinear case necessarily introduces overapproximation error. Consider the linear case first:

julia> using LazySets, TaylorModels

julia> const IA = IntervalArithmetic;

julia> I = IA.Interval(-0.5, 0.5) # interval remainder
[-0.5, 0.5]

julia> x₀ = IA.Interval(0.0) # expansion point
[0, 0]

julia> D = IA.Interval(-3.0, 1.0)
[-3, 1]

julia> p1 = Taylor1([2.0, 1.0], 2) # define a linear polynomial
 2.0 + 1.0 t + 𝒪(t³)

julia> p2 = Taylor1([0.9, 3.0], 2) # define another linear polynomial
 0.9 + 3.0 t + 𝒪(t³)

julia> vTM = [TaylorModel1(pi, I, x₀, D) for pi in [p1, p2]]
2-element Vector{TaylorModel1{Float64, Float64}}:
  2.0 + 1.0 t + [-0.5, 0.5]
  0.9 + 3.0 t + [-0.5, 0.5]

Here, vTM is a taylor model vector, since each component is a taylor model in one variable (TaylorModel1). Using overapproximate(vTM, Zonotope) we can compute its associated zonotope in generator representation:

julia> Z = overapproximate(vTM, Zonotope);

julia> center(Z)
2-element Vector{Float64}:
  1.0
 -2.1

julia> Matrix(genmat(Z))
2×3 Matrix{Float64}:
 2.0  0.5  0.0
 6.0  0.0  0.5

Note how the generators of this zonotope mainly consist of two pieces: one comes from the linear part of the polynomials, and another one that corresponds to the interval remainder. This conversion gives the same upper and lower bounds as the range evaluation using interval arithmetic:

julia> X = box_approximation(Z)
Hyperrectangle{Float64, Vector{Float64}, Vector{Float64}}([1.0, -2.1], [2.5, 6.5])

julia> Y = evaluate(vTM[1], vTM[1].dom) × evaluate(vTM[2], vTM[2].dom)
[-1.5, 3.5] × [-8.60001, 4.40001]

julia> H = convert(Hyperrectangle, Y) # this IntervalBox is the same as X
Hyperrectangle{Float64, StaticArraysCore.SVector{2, Float64}, StaticArraysCore.SVector{2, Float64}}([1.0, -2.1000000000000005], [2.5, 6.500000000000001])

However, the zonotope returns better results if we want to approximate the TM, since it is not axis-aligned:

julia> d = [-0.35, 0.93];

julia> ρ(d, Z) < ρ(d, X)
true

This function also works if the polynomials are non-linear; for example suppose that we add a third polynomial with a quadratic term:

julia> p3 = Taylor1([0.9, 3.0, 1.0], 3)
 0.9 + 3.0 t + 1.0 t² + 𝒪(t⁴)

julia> vTM = [TaylorModel1(pi, I, x₀, D) for pi in [p1, p2, p3]]
3-element Vector{TaylorModel1{Float64, Float64}}:
           2.0 + 1.0 t + [-0.5, 0.5]
           0.9 + 3.0 t + [-0.5, 0.5]
  0.9 + 3.0 t + 1.0 t² + [-0.5, 0.5]

julia> Z = overapproximate(vTM, Zonotope);

julia> center(Z)
3-element Vector{Float64}:
  1.0
 -2.1
  2.4

julia> Matrix(genmat(Z))
3×4 Matrix{Float64}:
 2.0  0.5  0.0  0.0
 6.0  0.0  0.5  0.0
 6.0  0.0  0.0  5.0

The fourth and last generator corresponds to the addition of the interval remainder and the box overapproximation of the nonlinear part of p3 over the domain.

Algorithm

Let $\text{vTM} = (p, I)$ be a vector of $m$ taylor models, where $I$ is the interval remainder in $\mathbb{R}^m$. Let $p_{lin}$ (resp. $p_{nonlin}$) correspond to the linear (resp. nonlinear) part of each scalar polynomial.

The range of $\text{vTM}$ can be enclosed by a zonotope with center $c$ and matrix of generators $G$, $Z = ⟨c, G⟩$, by performing a conservative linearization of $\text{vTM}$:

\[ vTM' = (p', I') := (p_{lin} − p_{nonlin} , I + \text{Int}(p_{nonlin})).\]

This algorithm proceeds in two steps:

1- Conservatively linearize $\text{vTM}$ as above and compute a box overapproximation of the nonlinear part. 2- Transform the linear taylor model to a zonotope exactly through variable normalization onto the symmetric intervals $[-1, 1]$.

source
overapproximate(vTM::Vector{TaylorModelN{N, T, S}}, ::Type{<:Zonotope};
                [remove_zero_generators]::Bool=true
                [normalize]::Bool=true) where {N, T, S}

Overapproximate a multivariate taylor model with a zonotope.

Input

  • vTMTaylorModelN
  • Zonotope – type for dispatch
  • remove_zero_generators – (optional; default: true) flag to remove zero generators of the resulting zonotope
  • normalize – (optional; default: true) flag to skip the normalization of the Taylor models

Output

A zonotope that overapproximates the range of the given taylor model.

Examples

Consider a vector of two 2-dimensional taylor models of order 2 and 4 respectively.

julia> using LazySets, TaylorModels

julia> const IA = IntervalArithmetic;

julia> x₁, x₂ = set_variables(Float64, ["x₁", "x₂"], order=8)
2-element Vector{TaylorN{Float64}}:
  1.0 x₁ + 𝒪(‖x‖⁹)
  1.0 x₂ + 𝒪(‖x‖⁹)

julia> x₀ = IntervalBox(0..0, 2) # expansion point
[0, 0]²

julia> Dx₁ = IA.Interval(0.0, 3.0) # domain for x₁
[0, 3]

julia> Dx₂ = IA.Interval(-1.0, 1.0) # domain for x₂
[-1, 1]

julia> D = Dx₁ × Dx₂ # take the Cartesian product of the domain on each variable
[0, 3] × [-1, 1]

julia> r = IA.Interval(-0.5, 0.5) # interval remainder
[-0.5, 0.5]

julia> p1 = 1 + x₁^2 - x₂
 1.0 - 1.0 x₂ + 1.0 x₁² + 𝒪(‖x‖⁹)

julia> p2 = x₂^3 + 3x₁^4 + x₁ + 1
 1.0 + 1.0 x₁ + 1.0 x₂³ + 3.0 x₁⁴ + 𝒪(‖x‖⁹)

julia> vTM = [TaylorModelN(pi, r, x₀, D) for pi in [p1, p2]]
2-element Vector{TaylorModelN{2, Float64, Float64}}:
            1.0 - 1.0 x₂ + 1.0 x₁² + [-0.5, 0.5]
  1.0 + 1.0 x₁ + 1.0 x₂³ + 3.0 x₁⁴ + [-0.5, 0.5]

julia> Z = overapproximate(vTM, Zonotope);

julia> center(Z)
2-element Vector{Float64}:
   5.5
 124.0

julia> Matrix(genmat(Z))
2×4 Matrix{Float64}:
 0.0  -1.0  5.0    0.0
 1.5   0.0  0.0  123.0

Algorithm

We refer to the algorithm description for the univariate case.

source
LazySets.Approximations._overapproximate_zonotope_vrepFunction
_overapproximate_zonotope_vrep(X::ConvexSet{N},
                               dir::AbstractDirections;
                               solver=default_lp_solver(N)) where {N}

Overapproximate a polytopic set with a zonotope of minimal total generator sum using only generators in the given directions.

Input

  • X – polytopic set
  • dir – directions used for the generators
  • solver – (optional, default: default_lp_solver(N)) the backend used to solve the linear program

Output

A zonotope that overapproximates X and uses at most the directions provided in dir (redundant directions will be ignored). Under this constraint, the zonotope has the minimal sum of generator vectors.

Notes

The algorithm only requires one representative of each generator direction and their additive inverse (e.g. only one of [1, 0] and [-1, 0]) and assumes that the directions are normalized. We preprocess the directions in that respect.

Algorithm

We solve a linear program parametric in the vertices $v_j$ of X and the directions $d_k$ in dir presented in Section 4.2 in [1], adapting the notation to the one used in this library.

\[ \min \sum_{k=1}^l α_k \ s.t. \ c + \sum_{k=1}^l b_{kj} * d_k = v_j \quad \forall j \ -α_k ≤ b_{kj} ≤ α_k \quad \forall k, j \ α_k ≥ 0 \quad \forall k\]

The resulting zonotope has center c and generators α_k · d_k.

Note that the first type of side constraints is vector-based and that the nonnegativity constraints (last type) are not stated explicitly in [1].

[1] Zonotopes as bounding volumes, L. J. Guibas et al, Proc. of Symposium on Discrete Algorithms, pp. 803-812.

source
LazySets.Approximations._overapproximate_zonotope_cpaFunction
_overapproximate_zonotope_cpa(X::ConvexSet, dir::Type{<:AbstractDirections})

Overapproximate a polytopic set with a zonotope using cartesian decomposition.

Input

  • X – polytopic set
  • dir – directions used for the generators

Output

A zonotope that overapproximates X.

Notes

The algorithm decomposes X in 2D sets and overapproximates those sets with zonotopes, and finally takes the cartesian product of the sets and converts to a zonotope.

Algorithm

The algorithm used is based on the section 8.2.4 of [1].

[1] Le Guernic, C. (2009). Reachability analysis of hybrid systems with linear continuous dynamics (Doctoral dissertation).

source

Underapproximations

LazySets.Approximations.underapproximateFunction
underapproximate(X::ConvexSet{N}, dirs::AbstractDirections;
                [apply_convex_hull]::Bool=false) where {N}

Compute the underapproximation of a convex set by sampling support vectors.

Input

  • X – set
  • dirs – directions
  • apply_convex_hull – (optional, default: false) if true, post-process the support vectors with a convex hull operation

Output

The VPolytope obtained by taking the convex hull of the support vectors of X along the directions determined by dirs.

Notes

Since the support vectors are not always unique, this algorithm may return a strict underapproximation even if the set can be exactly approximated using the given template.

source
underapproximate(X::ConvexSet, ::Type{<:Hyperrectangle};
                 solver=nothing) where {N}

Underapproximate a polygon with a hyperrectangle of maximal area.

Input

  • X – convex polygon
  • Hyperrectangle – type for dispatch
  • solver – (optional; default: nothing) nonlinear solver; if nothing, default_nln_solver(N) will be used

Output

A hyperrectangle underapproximation with maximal area.

Notes

The implementation only works for 2D sets, but the algorithm can be generalized.

Due to numerical issues, the result may be slightly outside the set.

Algorithm

The algorithm is taken from [1, Theorem 17] and solves a convex program (in fact a linear program with nonlinear objective). (The objective is modified to an equivalent version due to solver issues.)

[1] Mehdi Behroozi - Largest inscribed rectangles in geometric convex sets. arXiv:1905.13246.

source

Approximations

LazySets.Approximations.approximateFunction
approximate(R::Rectification; apply_convex_hull::Bool=false)

Approximate a rectification of a polytopic set with a convex polytope.

Input

  • R – rectification
  • apply_convex_hull – (optional; default: false) option to remove redundant vertices

Output

A polytope in vertex representation. There is no guarantee that the result over- or underapproximates R.

Algorithm

Let $X$ be the set that is rectified. We compute the vertices of $X$, rectify them, and return the convex hull of the result.

Notes

Let $X$ be the set that is rectified and let $p$ and $q$ be two vertices on a facet of $X$. Intuitively, an approximation may occur if the line segment connecting these vertices crosses a coordinate hyperplane and if the line segment connecting the rectified vertices has a different angle.

As a corollary, the approximation is exact for the special cases that the original set is contained in either the positive or negative orthant or is axis-aligned.

source

Box Approximations

LazySets.Approximations.ballinf_approximationFunction
ballinf_approximation(S::ConvexSet)

Overapproximate a set by a tight ball in the infinity norm.

Input

  • S – set

Output

A tight ball in the infinity norm.

Algorithm

The center and radius of the box are obtained by evaluating the support function of the given set along the canonical directions.

source
LazySets.Approximations.box_approximationFunction
box_approximation(S::ConvexSet{N}) where {N}

Overapproximate a set by a tight hyperrectangle.

Input

  • S – set

Output

A tight hyperrectangle.

Algorithm

The center of the hyperrectangle is obtained by averaging the support function of the given set in the canonical directions, and the lengths of the sides can be recovered from the distance among support functions in the same directions.

source
box_approximation(S::CartesianProductArray{N, <:AbstractHyperrectangle}) where {N}

Return a tight overapproximation of the Cartesian product array of a finite number of convex sets with and hyperrectangle.

Input

  • S – Cartesian product array of a finite number of convex set

Output

A hyperrectangle.

Algorithm

This method falls back to the corresponding convert method. Since the sets wrapped by the Cartesian product array are hyperrectangles, it can be done efficiently without overapproximation.

source
box_approximation(S::CartesianProduct{N, <:AbstractHyperrectangle, <:AbstractHyperrectangle}) where {N}

Return a tight overapproximation of the Cartesian product of two hyperrectangles by a new hyperrectangle.

Input

  • S – Cartesian product of two hyperrectangular sets

Output

A hyperrectangle.

Algorithm

This method falls back to the corresponding convert method. Since the sets wrapped by the Cartesian product are hyperrectangles, it can be done efficiently without overapproximation.

source
box_approximation(lm::LinearMap{N, <:AbstractHyperrectangle}) where {N}

Return a tight overapproximation of the linear map of a hyperrectangular set using a hyperrectangle.

Input

  • S – linear map of a hyperrectangular set

Output

A hyperrectangle.

Algorithm

If c and r denote the center and vector radius of a hyperrectangle H, a tight hyperrectangular overapproximation of M * H is obtained by transforming c ↦ M*c and r ↦ abs.(M) * r, where abs.(⋅) denotes the element-wise absolute value operator.

source
box_approximation(r::Rectification{N}) where {N}

Overapproximate the rectification of a convex set by a tight hyperrectangle.

Input

  • r – rectification of a convex set

Output

A hyperrectangle.

Algorithm

Box approximation and rectification distribute. Hence we first check whether the wrapped set is empty. If so, we return the empty set. Otherwise, we compute the box approximation of the wrapped set, rectify the resulting box (which is simple), and finally convert the resulting set to a box.

source
box_approximation(Z::AbstractZonotope)

Return a tight overapproximation of a zonotope with an axis-aligned box.

Input

  • Z – zonotope

Output

A hyperrectangle.

Algorithm

This function implements the method in [Section 5.1.2, 1]. A zonotope $Z = ⟨c, G⟩$ can be overapproximated tightly by an axis-aligned box (i.e. a Hyperrectangle) such that its center is $c$ and the radius along dimension $i$ is the column-sum of the absolute values of the $i$-th row of $G$ for $i = 1,…, p$, where $p$ is the number of generators of $Z$.

[1] Althoff, M., Stursberg, O., & Buss, M. (2010). Computing reachable sets of hybrid systems using a combination of zonotopes and polytopes. Nonlinear analysis: hybrid systems, 4(2), 233-249.

source
box_approximation(am::AbstractAffineMap{N, <:AbstractHyperrectangle}) where {N}

Overapproximate the affine map of a hyperrectangular set using a hyperrectangle.

Input

  • am – affine map of a hyperrectangular set

Output

A hyperrectangle.

Algorithm

If c and r denote the center and vector radius of a hyperrectangle H and v the translation vector, a tight hyperrectangular overapproximation of M * H + v is obtained by transforming c ↦ M*c+v and r ↦ abs.(M) * r, where abs.(⋅) denotes the element-wise absolute value operator.

source
box_approximation(ch::ConvexHull; [algorithm]::String="box")

Overapproximate a convex hull with a tight hyperrectangle.

Input

  • ch – convex hull
  • algorithm – (optional; default: "box") algorithm choice

Output

A hyperrectangle.

Algorithm

Let X and Y be the two sets of ch. We make use of the following property:

\[\square(CH(X, Y)) = \square\left( X \cup Y \right) = \square\left( \square(X) \cup \square(Y) \right)\]

If algorithm == "extrema", we compute the low and high coordinates of X and Y via extrema.

If algorithm == "box", we instead compute the box approximations of X and Y via box_approximation.

In both cases we then take the box approximation of the result.

The "extrema" algorithm is more efficient if extrema is efficient because it does not need to allocate the intermediate hyperrectangles.

source
box_approximation(ms::MinkowskiSum)

Compute the box approximation of the Minkowski sum of two sets.

Input

  • ms – Minkowski sum

Output

A hyperrectangle representing the box approximation of ms.

Algorithm

The box approximation distributes over the Minkowski sum:

\[\square(X \oplus Y) = \square(X) \oplus \square(Y)\]

It suffices to compute the box approximation of each summand and then take the concrete Minkowski sum for hyperrectangles.

source
LazySets.Approximations.symmetric_interval_hullFunction
symmetric_interval_hull(S::ConvexSet{N}) where {N}

Overapproximate a set by a tight hyperrectangle centered in the origin.

Input

  • S – set

Output

A tight hyperrectangle that is centrally symmetric wrt. the origin.

Algorithm

The center of the box is the origin, and the radius is obtained by computing the maximum value of the support function evaluated in the canonical directions.

Notes

The result is a hyperrectangle and hence in particular convex.

An alias for this function is box_approximation_symmetric.

source
LazySets.Approximations.box_approximation_helperFunction
box_approximation_helper(S::ConvexSet{N}) where {N}

Common code of box_approximation and box_approximation_symmetric.

Input

  • S – set

Output

A tuple containing the data that is needed to construct a tightly overapproximating hyperrectangle.

  • c – center
  • r – radius

Algorithm

The center of the hyperrectangle is obtained by averaging the support function of the given set in the canonical directions. The lengths of the sides can be recovered from the distance among support functions in the same directions.

source

Iterative refinement

LazySets.Approximations.LocalApproximationType
LocalApproximation{N, VN<:AbstractVector{N}}

Type that represents a local approximation in 2D.

Fields

  • p1 – first inner point
  • d1 – first direction
  • p2 – second inner point
  • d2 – second direction
  • q – intersection of the lines l1 ⟂ d1 at p1 and l2 ⟂ d2 at p2
  • refinable – states if this approximation is refinable
  • err – error upper bound

Notes

The criteria for being refinable are determined in the method new_approx.

source
LazySets.Approximations.PolygonalOverapproximationType
PolygonalOverapproximation{N, SN<:ConvexSet{N}, VN<:AbstractVector{N}}

Type that represents the polygonal approximation of a convex set.

Fields

  • S – convex set
  • approx_stack – stack of local approximations that still need to be examined
  • constraints – vector of linear constraints that are already finalized (i.e., they satisfy the given error bound)
source
LazySets.Approximations.new_approxMethod
new_approx(S::ConvexSet, p1::VN, d1::VN,
           p2::VN, d2::VN) where {N<:AbstractFloat, VN<:AbstractVector{N}}

Create a LocalApproximation instance for the given excerpt of a polygonal approximation.

Input

  • S – convex set
  • p1 – first inner point
  • d1 – first direction
  • p2 – second inner point
  • d2 – second direction

Output

A local approximation of S in the given directions.

source
LazySets.Approximations.addapproximation!Method
addapproximation!(Ω::PolygonalOverapproximation, p1::VN, d1::VN,
                  p2::VN, d2::VN) where {N, VN<:AbstractVector{N}}

Input

  • Ω – polygonal overapproximation of a convex set
  • p1 – first inner point
  • d1 – first direction
  • p2 – second inner point
  • d2 – second direction

Output

The list of local approximations in Ω of the set Ω.S is updated in-place and the new approximation is returned by this function.

source
LazySets.Approximations.refineMethod
refine(approx::LocalApproximation, S::ConvexSet)

Refine a given local approximation of the polygonal approximation of a convex set by splitting along the normal direction of the approximation.

Input

  • approx – local approximation to be refined
  • S – 2D convex set

Output

The tuple consisting of the refined right and left local approximations.

source
LazySets.Approximations.tohrepMethod
tohrep(Ω::PolygonalOverapproximation)

Convert a polygonal overapproximation into a concrete polygon.

Input

  • Ω – polygonal overapproximation of a convex set

Output

A polygon in constraint representation.

Algorithm

Internally we keep the constraints sorted. Hence we do not need to use addconstraint! when creating the HPolygon.

source
LazySets.Approximations._approximateMethod
_approximate(S::ConvexSet{N}, ε::Real) where {N<:AbstractFloat}

Return an ε-close approximation of the given 2D convex set (in terms of Hausdorff distance) as an inner and an outer approximation composed by sorted local Approximation2D.

Input

  • S – 2D convex set
  • ε – error bound

Output

An ε-close approximation of the given 2D convex set.

source
LazySets.Approximations.constraintMethod
constraint(approx::LocalApproximation)

Convert a local approximation to a linear constraint.

Input

  • approx – local approximation

Output

A linear constraint.

source

Template directions

LazySets.Approximations.AbstractDirectionsType
AbstractDirections{N, VN}

Abstract type for template direction representations.

Notes

This type is parameterzed by N and VN, where:

  • N stands for the numeric type
  • VN stands for the vector type with coefficients of type N

Each subtype is an iterator over a set of prescribed directions.

All subtypes should implement the standard iterator methods from Base, namely Base.length (returns the number of directions in the template), and Base.iterate. Moreover, the following methods should be implemented:

  • dim – return the ambient dimension of the template
  • eltype – return the type of each vector in the template

Optionally, subtypes may implement:

  • isbounding – (defaults to false) return true if an overapproximation with a list of template directions results in a bounded set, given a bounded input set, and false otherwise
  • isnormalized – (defaults to false) returns true if each direction in the given template has norm one w.r.t. the usual vector 2-norm
source
LazySets.Approximations.isboundingFunction
isbounding(ad::Type{<:AbstractDirections})

Checks if an overapproximation with a list of template directions results in a bounded set, given a bounded input set.

Input

  • ad – template directions

Output

Given a bounded set $X$, we can construct an outer approximation of $X$ by using the template directions ad as normal vectors of the facets. If this function returns true, then the result is again a bounded set (i.e., a polytope). Note that the result does not depend on the specific shape of $X$, as long as $X$ is bounded.

Notes

By default, this function returns false in order to be conservative. Custom subtypes of AbstractDirections should hence add a method for this function.

source
LazySets.Approximations.isnormalizedFunction
isnormalized(ad::Type{<:AbstractDirections})

Returns whether the given template directions is normalized with respect to the 2-norm.

Input

  • ad – template directions

Output

true if the 2-norm of each element in ad is one and false otherwise.

source
LazySets.projectMethod
project(S::ConvexSet,
        block::AbstractVector{Int},
        directions::Type{<:AbstractDirections},
        [n]::Int;
        [kwargs...]
       )

Project a high-dimensional set to a given block using template directions.

Input

  • S – set
  • block – block structure - a vector with the dimensions of interest
  • directions – template directions
  • n – (optional, default: dim(S)) ambient dimension of the set S

Output

The template direction approximation of the projection of S.

source
LazySets.Approximations.BoxDirectionsType
BoxDirections{N, VN} <: AbstractDirections{N, VN}

Box directions representation.

Fields

  • n – dimension

Notes

Box directions can be seen as the vectors where only one entry is ±1, and all other entries are 0. In dimension $n$, there are $2n$ such directions.

The default vector representation used in this template is a ReachabilityBase.Arrays.SingleEntryVector, although other implementations can be used such as a regular Vector and a sparse vector, SparseVector.

Examples

The template can be constructed by passing the dimension. For example, in dimension two,

julia> dirs = BoxDirections(2)
BoxDirections{Float64, ReachabilityBase.Arrays.SingleEntryVector{Float64}}(2)

julia> length(dirs)
4

By default, each direction is represented in this iterator as a SingleEntryVector, i.e. a vector with only one non-zero element,

julia> eltype(dirs)
ReachabilityBase.Arrays.SingleEntryVector{Float64}

In two dimensions, the directions defined by BoxDirections are normal to the facets of a box.

julia> collect(dirs)
4-element Vector{ReachabilityBase.Arrays.SingleEntryVector{Float64}}:
 [1.0, 0.0]
 [0.0, 1.0]
 [0.0, -1.0]
 [-1.0, 0.0]

The numeric type can be specified as well:

julia> BoxDirections{Rational{Int}}(10)
BoxDirections{Rational{Int64}, ReachabilityBase.Arrays.SingleEntryVector{Rational{Int64}}}(10)

julia> length(ans)
20
source
LazySets.Approximations.DiagDirectionsType
DiagDirections{N, VN} <: AbstractDirections{N, VN}

Diagonal directions representation.

Fields

  • n – dimension

Notes

Diagonal directions can be seen as all diagonal directions (all entries are ±1). In dimension $n$, there are in total $2^n$ such directions.

Examples

The template can be constructed by passing the dimension. For example, in dimension two,

julia> dirs = DiagDirections(2)
DiagDirections{Float64, Vector{Float64}}(2)

julia> length(dirs) # number of directions
4

By default, each direction is represented in this iterator as a regular vector:

julia> eltype(dirs)
Vector{Float64} (alias for Array{Float64, 1})

In two dimensions, the directions defined by DiagDirections are normal to the facets of a ball in the 1-norm.

julia> collect(dirs)
4-element Vector{Vector{Float64}}:
 [1.0, 1.0]
 [-1.0, 1.0]
 [1.0, -1.0]
 [-1.0, -1.0]

The numeric type can be specified as well:

julia> DiagDirections{Rational{Int}}(10)
DiagDirections{Rational{Int64}, Vector{Rational{Int64}}}(10)

julia> length(ans)
1024
source
LazySets.Approximations.OctDirectionsType
OctDirections{N, VN} <: AbstractDirections{N, VN}

Octagon directions representation.

Fields

  • n – dimension

Notes

Octagon directions consist of all vectors that are zero almost everywhere except in two dimensions $i$, $j$ (possibly $i = j$) where it is $±1$. In dimension $n$, there are $2n^2$ such directions.

Examples

The template can be constructed by passing the dimension. For example, in dimension two,

julia> dirs = OctDirections(2)
OctDirections{Float64, SparseArrays.SparseVector{Float64, Int64}}(2)

julia> length(dirs) # number of directions
8

By default, each direction is represented in this iterator as a sparse vector:

julia> eltype(dirs)
SparseArrays.SparseVector{Float64, Int64}

In two dimensions, the directions defined by OctDirections are normal to the facets of an octagon.

julia> first(dirs)
2-element SparseArrays.SparseVector{Float64, Int64} with 2 stored entries:
  [1]  =  1.0
  [2]  =  1.0

julia> Vector.(collect(dirs))
8-element Vector{Vector{Float64}}:
 [1.0, 1.0]
 [1.0, -1.0]
 [-1.0, 1.0]
 [-1.0, -1.0]
 [1.0, 0.0]
 [0.0, 1.0]
 [0.0, -1.0]
 [-1.0, 0.0]

The numeric type can be specified as well:

julia> OctDirections{Rational{Int}}(10)
OctDirections{Rational{Int64}, SparseArrays.SparseVector{Rational{Int64}, Int64}}(10)

julia> length(ans)
200
source
LazySets.Approximations.BoxDiagDirectionsType
BoxDiagDirections{N, VN} <: AbstractDirections{N, VN}

Box-diagonal directions representation.

Fields

  • n – dimension

Notes

Box-diagonal directions can be seen as the union of diagonal directions (all entries are ±1) and box directions (one entry is ±1, all other entries are 0). The iterator first enumerates all diagonal directions, and then all box directions. In dimension $n$, there are in total $2^n + 2n$ such directions.

Examples

The template can be constructed by passing the dimension. For example, in dimension two,

julia> dirs = BoxDiagDirections(2)
BoxDiagDirections{Float64, Vector{Float64}}(2)

julia> length(dirs) # number of directions
8

By default, each direction is represented in this iterator as a regular vector:

julia> eltype(dirs)
Vector{Float64} (alias for Array{Float64, 1})

In two dimensions, the directions defined by BoxDiagDirections are normal to the facets of an octagon.

julia> collect(dirs)
8-element Vector{Vector{Float64}}:
 [1.0, 1.0]
 [-1.0, 1.0]
 [1.0, -1.0]
 [-1.0, -1.0]
 [1.0, 0.0]
 [0.0, 1.0]
 [0.0, -1.0]
 [-1.0, 0.0]

The numeric type can be specified as well:

julia> BoxDiagDirections{Rational{Int}}(10)
BoxDiagDirections{Rational{Int64}, Vector{Rational{Int64}}}(10)

julia> length(ans)
1044
source
LazySets.Approximations.PolarDirectionsType
PolarDirections{N<:AbstractFloat, VN<:AbstractVector{N}} <: AbstractDirections{N, VN}

Polar directions representation.

Fields

  • – length of the partition of the polar angle

Notes

The PolarDirections constructor provides a sample of the unit sphere in $\mathbb{R}^2$, which is parameterized by the polar angle $φ ∈ Dφ := [0, 2π]$; see the wikipedia entry Polar coordinate system for details.

The integer argument $Nφ$ defines how many samples of $Dφ$ are taken. The Cartesian components of each direction are obtained with

\[[cos(φᵢ), sin(φᵢ)].\]

Examples

The integer passed as an argument is used to discretize $φ$:

julia> pd = PolarDirections(2);

julia> pd.stack
2-element Vector{Vector{Float64}}:
 [1.0, 0.0]
 [-1.0, 1.2246467991473532e-16]

julia> length(pd)
2
source
LazySets.Approximations.SphericalDirectionsType
SphericalDirections{N<:AbstractFloat, VN<:AbstractVector{N}} <: AbstractDirections{N, VN}

Spherical directions representation.

Fields

  • – length of the partition of the azimuthal angle
  • – length of the partition of the polar angle
  • stack – list of computed directions

Notes

The SphericalDirections constructor provides a sample of the unit sphere in $\mathbb{R}^3$, which is parameterized by the azimuthal and polar angles $θ ∈ Dθ := [0, π]$ and $φ ∈ Dφ := [0, 2π]$ respectively, see the wikipedia entry Spherical coordinate system for details.

The integer arguments $Nθ$ and $Nφ$ define how many samples along the domains $Dθ$ and $Dφ$ respectively are taken. The Cartesian components of each direction are obtained with

\[[sin(θᵢ)*cos(φᵢ), sin(θᵢ)*sin(φᵢ), cos(θᵢ)].\]

The north and south poles are treated separately so that those points are not considered more than once.

Examples

A SphericalDirections template can be built in different ways. If you pass only one integer, the same value is used to discretize both $θ$ and $φ$:

julia> sd = SphericalDirections(3);

julia> sd.Nθ, sd.Nφ
(3, 3)

julia> length(sd)
4

Pass two integers to control the discretization in $θ$ and in $φ$ separately:

julia> sd = SphericalDirections(4, 5);

julia> length(sd)
10

julia> sd = SphericalDirections(4, 8);

julia> length(sd)
16
source
LazySets.Approximations.CustomDirectionsType
CustomDirections{N, VN<:AbstractVector{N}} <: AbstractDirections{N, VN}

User-defined template directions.

Fields

  • directions – list of template directions
  • n – (optional; default: computed from `directions) dimension
  • check_boundedness – (optional; default: true) flag to check boundedness
  • check_normalization – (optional; default: true) flag to check whether all directions are normalized

Notes

This struct is a wrapper type for a set of user-defined directions which are iterated over. It has fields for the list of directions, the set dimension, and (boolean) cache fields for the boundedness and normalization properties. The latter are checked by default upon construction.

To check boundedness, we overapproximate the unit ball in the infinity norm using the given directions and check if the resulting set is bounded.

The dimension will also be determined automatically, unless the empty vector is passed (in which case the optional argument n needs to be specified).

Examples

Creating a template with box directions in dimension two:

julia> dirs = CustomDirections([[1.0, 0.0], [-1.0, 0.0], [0.0, 1.0], [0.0, -1.0]]);

julia> dirs.directions
4-element Vector{Vector{Float64}}:
 [1.0, 0.0]
 [-1.0, 0.0]
 [0.0, 1.0]
 [0.0, -1.0]

julia> LazySets.Approximations.isbounding(dirs)
true

julia> LazySets.Approximations.isnormalized(dirs)
true
source

See also overapproximate(X::ConvexSet, dir::AbstractDirections)::HPolytope.

Hausdorff distance

LazySets.Approximations.hausdorff_distanceFunction
hausdorff_distance(X::ConvexSet{N}, Y::ConvexSet{N}; [p]::N=N(Inf),
                   [ε]=N(1e-3)) where {N}

Compute the Hausdorff distance between two convex sets up to a given threshold.

Input

  • X – convex set
  • Y – convex set
  • p – (optional, default: Inf) norm parameter of the Hausdorff distance
  • ε – (optional, default: 1e-3) precision threshold; the true Hausdorff distance may diverge from the result by at most this value

Output

A value from the $ε$-neighborhood of the Hausdorff distance between $X$ and $Y$.

Notes

Given a $p$-norm, the Hausdorff distance $d_H^p(X, Y)$ between sets $X$ and $Y$ is defined as follows:

\[ d_H^p(X, Y) = \inf\{δ ≥ 0 \mid Y ⊆ X ⊕ δ 𝐵_p^n \text{ and } X ⊆ Y ⊕ δ 𝐵_p^n\}\]

Here $𝐵_p^n$ is the $n$-dimensional unit ball in the $p$-norm.

The implementation may internally rely on the support function of $X$ and $Y$; hence any imprecision in the implementation of the support function may affect the result. At the time of writing, the only set type with imprecise support function is the lazy Intersection.

Algorithm

We perform binary search for bounding the Hausdorff distance in an interval $[l, u]$, where initially $l$ is $0$ and $u$ is described below. The binary search terminates when $u - l ≤ ε$, i.e., the interval becomes sufficiently small.

To find an upper bound $u$, we start with the heuristics of taking the biggest distance in the axis-parallel directions. As long as this bound does not work, we increase the bound by $2$.

Given a value $δ$, to check whether the sets are within Hausdorff distance $δ$, we simply check the inclusions given above, where on the right-hand side we use a lazy Bloating.

source