Utility functions

Module Arrays.jl – Auxiliary machinery for vectors and matrices.

cross_product(M::AbstractMatrix{N}) where {N<:Real}

Compute the high-dimensional cross product of $n-1$ $n$-dimensional vectors.

Input

• M – $n × n - 1$-dimensional matrix

Output

A vector.

Algorithm

The cross product is defined as follows:

\[\left[ \dots, (-1)^{n+1} \det(M^{[i]}), \dots \right]^T\]

where $M^{[i]}$ is defined as $M$ with the $i$-th row removed. See Althoff, Stursberg, Buss: Computing Reachable Sets of Hybrid Systems Using a Combination of Zonotopes and Polytopes. 2009.

nonzero_columns(A::AbstractMatrix)

Return all columns that have at least one non-zero entry.

Input

• A – matrix

Output

A vector of indices.

dot_zero(x::AbstractVector{N}, y::AbstractVector{N}) where{N<:Real}

Dot product with preference for zero value in the presence of infinity values.

Input

• x – first vector
• y – second vector

Output

The dot product of x and y, but with the rule that 0 * Inf == 0.

hasfullrowrank(M::AbstractMatrix)

Check whether a matrix has full row rank.

Input

• M – matrix

Output

true iff the matrix has full row rank.

inner(x::AbstractVector{N}, A::AbstractMatrix{N}, y::AbstractVector{N}
     ) where {N}

Compute the inner product $xᵀ A y$.

Input

• x – vector on the left
• A – matrix
• y – vector on the right

Output

The (scalar) result of the multiplication.

is_cyclic_permutation(candidate::AbstractVector, paragon::AbstractVector)

Checks if the elements in candidate are a cyclic permutation of the elements in paragon.

Input

• candidate – candidate vector
• paragon – paragon vector

Output

A boolean indicating if the elements of candidate are in the same order as in paragon or any of its cyclic permutations.

isinvertible(M::Matrix; [cond_tol]::Number=DEFAULT_COND_TOL)

A sufficient check of a matrix being invertible (or nonsingular).

Input

• M – matrix
• cond_tol – (optional, default: DEFAULT_COND_TOL) tolerance of matrix condition

Output

If the result is true, M is invertible. If the result is false, the matrix is non-square or this function could not conclude.

Algorithm

We check whether the matrix is square and whether the matrix condition number cond(M) is below some prescribed tolerance.

ismultiple(u::AbstractVector{<:Real}, v::AbstractVector{<:Real})

Check whether two vectors are linearly dependent.

Input

• u – first vector
• v – second vector

Output

(true, k) iff the vectors are identical up to a scaling factor k ≠ 0 such that u = k * v, and (false, 0) otherwise.

Examples

julia> using LazySets: ismultiple

julia> ismultiple([1, 2, 3], [2, 4, 6])
(true, 0.5)

julia> ismultiple([1, 2, 3], [3, 2, 1])
(false, 0)

julia> ismultiple([1, 2, 3], [-1, -2, -3])
(true, -1.0)

ispermutation(u::AbstractVector{T}, v::AbstractVector) where {T}

Check that two vectors contain the same elements up to reordering.

Input

• u – first vector
• v – second vector

Output

true iff the vectors are identical up to reordering.

Examples

julia> LazySets.ispermutation([1, 2, 2], [2, 2, 1])
true

julia> LazySets.ispermutation([1, 2, 2], [1, 1, 2])
false

Notes

Containment check is performed using LazySets._in(e, v), so in the case of floating point numbers, the precision to which the check is made is determined by the type of elements in v. See _in and _isapprox for more information.

Note that approximate equality is not an equivalence relation. Hence the result may depend on the order of the elements.

is_right_turn([O::AbstractVector{N}=[0, 0]], u::AbstractVector{N},
              v::AbstractVector{N}) where {N<:Real}

Determine whether the acute angle defined by three 2D points O, u, v in the plane is a right turn (< 180° counter-clockwise) with respect to the center O. Determine if the acute angle defined by two 2D vectors is a right turn (< 180° counter-clockwise) with respect to the center O.

Input

• O – (optional; default: [0, 0]) 2D center point
• u – first 2D direction
• v – second 2D direction

Output

true iff the vectors constitute a right turn.

issquare(M::AbstractMatrix)

Check whether a matrix is square.

Input

• M – matrix

Output

true iff the matrix is square.

nonzero_indices(v::AbstractVector{N}) where {N<:Real}

Return the indices in which a vector is non-zero.

Input

• v – vector

Output

A vector of ascending indices i such that the vector is non-zero in dimension i.

remove_duplicates_sorted!(v::AbstractVector)

Remove duplicate entries in a sorted vector.

Input

• v – sorted vector

Output

The input vector without duplicates.

remove_zero_columns(A::AbstractMatrix)

Return a matrix with all columns containing only zero entries removed.

Input

• A – matrix

Output

The original matrix A if it contains no zero columns or otherwise a new matrix where those columns have been removed.

right_turn([O::AbstractVector{N}=[0, 0]], u::AbstractVector{N},
           v::AbstractVector{N}) where {N<:Real}

Compute a scalar that determines whether the acute angle defined by three 2D points O, u, v in the plane is a right turn (< 180° counter-clockwise) with respect to the center O.

Input

• O – (optional; default: [0, 0]) 2D center point
• u – first 2D point
• v – second 2D point

Output

A scalar representing the rotation. If the result is 0, the points are collinear; if it is positive, the points constitute a positive angle of rotation around O from u to v; otherwise they constitute a negative angle.

Algorithm

The cross product is used to determine the sense of rotation.

samedir(u::AbstractVector{<:Real}, v::AbstractVector{<:Real})

Check whether two vectors point in the same direction.

Input

• u – first vector
• v – second vector

Output

(true, k) iff the vectors are identical up to a positive scaling factor k such that u = k * v, and (false, 0) otherwise.

Examples

julia> using LazySets: samedir

julia> samedir([1, 2, 3], [2, 4, 6])
(true, 0.5)

julia> samedir([1, 2, 3], [3, 2, 1])
(false, 0)

julia> samedir([1, 2, 3], [-1, -2, -3])
(false, 0)

SingleEntryVector{N} <: AbstractVector{N}

A sparse unit vector with arbitrary one-element.

Fields

• i – index of non-zero entry
• n – vector length
• v – non-zero entry

to_negative_vector(v::AbstractVector{N}) where {N}

Negate a vector and convert to type Vector.

Input

• v – vector

Output

A Vector equivalent to $-v$.

_up(u::AbstractVector, v::AbstractVector)

Checks if the given vector is pointing towards the given direction.

Input

• u – direction
• v – vector

Output

A boolean indicating if the vector is pointing towards the direction.

_dr(u::AbstractVector, Vi::AbstractVector, Vj::AbstractVector)

Returns the direction of the difference of the given vectors.

Input

• u – direction
• Vi – first vector
• Vj – second vector

Output

A number indicating the direction of the difference of the given vectors.

_above(u::AbstractVector, Vi::AbstractVector, Vj::AbstractVector)

Checks if the difference of the given vectors is pointing towards the given direction.

Input

• u – direction
• Vi – first vector
• Vj – second vector

Output

A boolean indicating if the difference of the given vectors is pointing towards the given direction.

minmax(A, B, C)

Compute the minimum and maximum of three numbers A, B, C.

Input

• A – first number
• B – second number
• C – third number

Output

The minimum and maximum of three given numbers.

Examples

julia> LazySets.minmax(1.4, 52.4, -5.2)
(-5.2, 52.4)

arg_minmax(A, B, C)

Compute the index (1, 2, 3) of the minimum and maximum of three numbers A, B, C.

Input

• A – first number
• B – second number
• C – third number

Output

The index of the minimum and maximum of the three given numbers.

Examples

julia> LazySets.arg_minmax(1.4, 52.4, -5.2)
(3, 2)

extend(M::AbstractMatrix; check_rank=true)

Return an invertible extension of M whose first n columns span the column space of M, assuming that size(M) = (m, n), m > n and the rank of M is n.

Input

• M – rectangular m × n matrix with m > n and full rank (i.e. its rank is n)
• check_rank – (optional, default: true) if true, check the rank assumption, otherwise do not perform this check

Output

The tuple (Mext, inv_Mext), where Mext is a square m × m invertible matrix that extends M, i.e. in the sense that Mext = [M | Q2], and the rank of Mext is m. Here, inv_Mext is the inverse of Mext.

Algorithm

First we compute the QR decomposition of M to extract a suitable subspace of column vectors (Q2) that are orthogonal to the column span of M. Then we observe that the inverse of the extended matrix Mext = [M | Q2] is [R⁻¹Qᵀ; Q2ᵀ].

projection_matrix(block::AbstractVector{Int}, n::Int, [N]::Type{<:Number}=Float64)

Return the projection matrix associated to the given block of variables.

Input

• block – integer vector with the variables of interest
• n – integer representing the ambient dimension
• N – (optional, default: Float64) number type

Output

A sparse matrix that corresponds to the projection onto the variables in block.

Examples

julia> using LazySets: projection_matrix

julia> projection_matrix([1, 3], 4)
2×4 SparseArrays.SparseMatrixCSC{Float64,Int64} with 2 stored entries:
  [1, 1]  =  1.0
  [2, 3]  =  1.0

julia> Matrix(ans)
2×4 Array{Float64,2}:
 1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0

_vector_type(T)

Return a corresponding vector type with respect to type T.

Input

• T – vector or matrix type

Output

A vector type that corresponds in some sense (see Notes below) to T.

Notes

• If T is a sparse vector or a sparse matrix, the corresponding type is also a sparse vector.
• If T is a regular vector (i.e. Vector) or a regular matrix (i.e. Matrix), the corresponding type is also a regular vector.
• Otherwise, the corresponding type is a regular vector.

_matrix_type(T)

Return a corresponding matrix type with respect to type T.

Input

• T – vector or matrix type

Output

A matrix type that corresponds in some sense (see Notes below) to T.

Notes

• If T is a sparse vector or a sparse matrix, the corresponding type is also a sparse matrix.
• If T is a regular vector (i.e. Vector) or a regular matrix (i.e. Matrix), the corresponding type is also a regular matrix.
• Otherwise, the corresponding type is a regular matrix.

allequal(x)

Check whether all elements in a sequence are equal

Input

• x – sequence

Output

true iff all elements in x are equal.

Notes

The code is taken from here.

distance(x::AbstractVector, y::AbstractVector, p::Real=2.0)

Compute the distance between two vectors with respect to the given p-norm, computed as

\[ \|x - y\|_p = \left( \sum_{i=1}^n | x_i - y_i |^p \right)^{1/p}\]

Input

• x – vector
• y – vector
• p – (optional, default: 2.0) the p-norm used; p = 2.0 corresponds to the usual Euclidean norm

Output

A scalar representing $‖ x - y ‖_p$.

same_sign(A::AbstractArray{N}; [optimistic]::Bool=false) where {N}

Check whether all elements of the given array have the same sign.

Input

• A – array
• optimistic – (optional; default: false) flag for expressing that the expected result is true

Output

true if and only if all elements in M have the same sign.

Algorithm

If optimistic is false, we check the sign of the first element and compare to the sign of all elements.

If optimistic is true, we compare the absolute element sum with the sum of the absolute of the elements; this is faster if the result is true because there is no branching.

\[ |\sum_i A_i| = \sum_i |A_i|\]

to_matrix(vectors::AbstractVector{VN},
          [m]=length(vectors[1]),
          [MT]=_matrix_type(VN)) where {VN}

Input

• vectors – list of vectors
• m – (optional; default: length(vectors[1])) number of rows
• MT – (optional; default: _matrix_type(VN)) type of target matrix

Output

A matrix with the column vectors from vectors in the same order.

_rationalize(::Type{T}, x::AbstractVecOrMat{N}, tol::Real) where {T<:Integer, N<:AbstractFloat}

Approximate an array of floating point numbers as a rational vector with entries of the given integer type.

Input

• T – (optional, default: Int) integer type to represent the rationals
• x – vector of floating point numbers
• tol – (optional, default: eps(N)) tolerance the result at entry i will differ from x[i] by no more than tol

Output

An array of type Rational{T} where the i-th entry is the rationalization of the i-th component of x.

Notes

See also Base.rationalize.

an_element_helper(hp::Hyperplane{N},
                  [nonzero_entry_a]::Int) where {N}

Helper function that computes an element on a hyperplane's hyperplane.

Input

• hp – hyperplane
• nonzero_entry_a – (optional, default: computes the first index) index i such that hp.a[i] is different from 0

Output

An element on a hyperplane.

Algorithm

We compute the point on the hyperplane as follows:

• We already found a nonzero entry of $a$ in dimension, say, $i$.
• We set $x[i] = b / a[i]$.
• We set $x[j] = 0$ for all $j ≠ i$.

binary_search_constraints(d::AbstractVector{N},
                          constraints::Vector{<:LinearConstraint{N}},
                          n::Int,
                          k::Int;
                          [choose_lower]::Bool=false) where {N}

Performs a binary search in the constraints.

Input

• d – direction
• constraints – constraints
• n – number of constraints
• k – start index
• choose_lower – (optional, default: false) flag for choosing the lower index (see the 'Output' section)

Output

In the default setting, the result is the smallest index k such that d <= constraints[k], or n+1 if no such k exists. If the choose_lower flag is set, the result is the largest index k such that constraints[k] < d, which is equivalent to being k-1 in the normal setting.

get_radius!(sih::SymmetricIntervalHull{N},
            i::Int,
            n::Int=dim(sih)) where {N}

Compute the radius of the symmetric interval hull of a set in a given dimension.

Input

• sih – symmetric interval hull of a set
• i – dimension in which the radius should be computed
• n – (optional, default: dim(sih)) set dimension

Output

The radius of the symmetric interval hull of a set in a given dimension.

Algorithm

We ask for the support vector of the underlying set for both the positive and negative unit vector in the dimension i.

is_tighter_same_dir_2D(c1::LinearConstraint,
                       c2::LinearConstraint;
                       [strict]::Bool=false)

Check if the first of two two-dimensional constraints with equivalent normal direction is tighter.

Input

• c1 – first linear constraint
• c2 – second linear constraint
• strict – (optional; default: false) check for strictly tighter constraints?

Output

true iff the first constraint is tighter.

sign_cadlag(x::N) where {N<:Real}

This function works like the sign function but is $1$ for input $0$.

Input

• x – real scalar

Output

$1$ if $x ≥ 0$, $-1$ otherwise.

Notes

This is the sign function right-continuous at zero (see càdlàg function). It can be used with vector-valued arguments via the dot operator.

Examples

julia> LazySets.sign_cadlag.([-0.6, 1.3, 0.0])
3-element Array{Float64,1}:
 -1.0
  1.0
  1.0

_leq_trig(u::AbstractVector{N}, v::AbstractVector{N}) 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.

_random_zero_sum_vector(rng::AbstractRNG, N::Type{<:Real}, n::Int)

Create a random vector with entries whose sum is zero.

Input

• rng – random number generator
• N – numeric type
• n – length of vector

Output

A random vector of random numbers such that all positive entries come first and all negative entries come last, and such that the total sum is zero.

Algorithm

This is a preprocessing step of the algorithm here based on P. Valtr. Probability that n random points are in convex position.

rectify(x::AbstractVector{N}) where {N<:Real}

Rectify a vector, i.e., take the element-wise maximum with zero.

Input

• x – vector

Output

A copy of the vector where each negative entry is replaced by zero.

rectify(X::LazySet, [concrete_intersection]::Bool=false)

Concrete rectification of a set.

Input

• X – set
• concrete_intersection – (optional, default: false) flag to compute concrete intersections for intermediate results

Output

A set corresponding to the rectification of X, which is in general a union of linear maps of intersections.

Algorithm

For each dimension in which X is both positive and negative we split X into these two parts. Additionally we project the negative part to zero.

rectify(H::AbstractHyperrectangle)

Concrete rectification of a hyperrectangular set.

Input

• H – hyperrectangular set

Output

The Hyperrectangle that corresponds to the rectification of H.

rectify(x::Interval{N}) where {N}

Concrete rectification of an interval.

Input

• x – interval

Output

The Interval that corresponds to the rectification of x.

Notes

Note that the result is an Interval even if the set becomes a singleton (which is the case if the original interval was nonpositive).

rectify(S::Singleton)

Concrete rectification of a singleton.

Input

• S – singleton

Output

The Singleton that corresponds to the rectification of S.

rectify(Z::ZeroSet)

Concrete rectification of a zero set.

Input

• Z – zero set

Output

The same set.

require(package::Symbol; fun_name::String="", explanation::String="")

Helper method to check for optional packages and printing an error message.

Input

• package – symbol of the package name
• fun_name – (optional; default: "") name of the function that requires the package
• explanation – (optional; default: "") additional explanation in the error message

Output

If the package is loaded, this function has no effect. Otherwise it prints an error message.

Algorithm

This function uses @assert and hence loses its ability to print an error message if assertions are deactivated.

reseed(rng::AbstractRNG, seed::Union{Int, Nothing})

Reset the RNG seed if the seed argument is a number.

Input

• rng – random number generator
• seed – seed for reseeding

Output

The input RNG if the seed is nothing, and a reseeded RNG otherwise.

sameblockstructure(x::AbstractVector{S1}, y::AbstractVector{S2} ) where {S1<:LazySet, S2<:LazySet}

Check whether two vectors of sets have the same block structure, i.e., the $i$-th entry in the vectors have the same dimension.

Input

• x – first vector
• y – second vector

Output

true iff the vectors have the same block structure.

substitute(substitution::Dict{Int, T}, x::AbstractVector{T}) where {T}

Apply a substitution to a given vector.

Input

• substitution – substitution (a mapping from an index to a new value)
• x – vector

Output

A fresh vector corresponding to x after substitution was applied.

substitute!(substitution::Dict{Int, T}, x::AbstractVector{T}) where {T}

Apply a substitution to a given vector.

Input

• substitution – substitution (a mapping from an index to a new value)
• x – vector (modified in this function)

Output

The same (but see the Notes below) vector x but after substitution was applied.

Notes

The vector x is modified in-place if it has type Vector or SparseVector. Otherwise, we first create a new Vector from it.

    σ_helper(d::AbstractVector,
             hp::Hyperplane;
             error_unbounded::Bool=true,
             [halfspace]::Bool=false)

Return the support vector of a hyperplane.

Input

• d – direction
• hp – hyperplane
• error_unbounded – (optional, default: true) true if an error should be thrown whenever the set is unbounded in the given direction
• halfspace – (optional, default: false) true if the support vector should be computed for a half-space

Output

A pair (v, b) where v is a vector and b is a Boolean flag.

The flag b is false in one of the following cases:

1. The direction has norm zero.
2. The direction is the hyperplane's normal direction.
3. The direction is the opposite of the hyperplane's normal direction and

halfspace is false. In all these cases, v is any point on the hyperplane.

Otherwise, the flag b is true, the set is unbounded in the given direction, and v is any vector.

If error_unbounded is true and the set is unbounded in the given direction, this function throws an error instead of returning.

Notes

For correctness, consider the weak duality of LPs: If the primal is unbounded, then the dual is infeasible. Since there is only a single constraint, the feasible set of the dual problem is hp.a ⋅ y == d, y >= 0 (with objective function hp.b ⋅ y). It is easy to see that this problem is infeasible whenever a is not parallel to d.

get_constrained_lowdimset(cpa::CartesianProductArray{N, S},
                          P::AbstractPolyhedron{N}
                          ) where {N, S<:LazySet{N}}

Preprocessing step for the intersection between a Cartesian product array and a polyhedron.

Input

• cpa – Cartesian product array of sets
• P – polyhedron

Output

A tuple of low-dimensional set, list of constrained dimensions, original block structure of low-dimensional set and corresponding blocks indices.

Notes

This method returns a four-tuple with the following entries:

1. a low-dimensional CartesianProductArray in the constrained dimensions of

the original cpa

1. the constrained variables and variables in corresponding blocks
2. the original block structure of the low-dimensional sets
3. the list of the constrained blocks.

@neutral(SET, NEUT)

Create functions to make a lazy set operation commutative with a given neutral element set type.

Input

• SET – lazy set operation type
• NEUT – set type for neutral element

Output

Nothing.

Notes

This macro generates four functions (possibly two more if @absorbing has been used in advance) (possibly two or four more if @declare_array_version has been used in advance).

Examples

@neutral(MinkowskiSum, N) creates at least the following functions:

• neutral(::MinkowskiSum) = N
• MinkowskiSum(X, N) = X
• MinkowskiSum(N, X) = X
• MinkowskiSum(N, N) = N

@absorbing(SET, ABS)

Create functions to make a lazy set operation commutative with a given absorbing element set type.

Input

• SET – lazy set operation type
• ABS – set type for absorbing element

Output

Nothing.

Notes

This macro generates four functions (possibly two more if @neutral has been used in advance) (possibly two or four more if @declare_array_version has been used in advance).

Examples

@absorbing(MinkowskiSum, A) creates at least the following functions:

• absorbing(::MinkowskiSum) = A
• MinkowskiSum(X, A) = A
• MinkowskiSum(A, X) = A
• MinkowskiSum(A, A) = A

@neutral_absorbing(SET, NEUT, ABS)

Create two functions to avoid method ambiguities for a lazy set operation with respect to neutral and absorbing element set types.

Input

• SET – lazy set operation type
• NEUT – set type for neutral element
• ABS – set type for absorbing element

Output

A quoted expression containing the function definitions.

Examples

@neutral_absorbing(MinkowskiSum, N, A) creates the following functions as quoted expressions:

• MinkowskiSum(N, A) = A
• MinkowskiSum(A, N) = A

@declare_array_version(SET, SETARR)

Create functions to connect a lazy set operation with its array set type.

Input

• SET – lazy set operation type
• SETARR – array set type

Output

Nothing.

Notes

This macro generates six functions (and possibly up to eight more if @neutral/@absorbing has been used in advance for the base and/or array set type).

Examples

@declare_array_version(MinkowskiSum, MinkowskiSumArray) creates at least the following functions:

• array_constructor(::MinkowskiSum) = MinkowskiSumArray
• is_array_constructor(::MinkowskiSumArray) = true
• MinkowskiSum!(X, Y)
• MinkowskiSum!(X, arr)
• MinkowskiSum!(arr, X)
• MinkowskiSum!(arr1, arr2)

@array_neutral(FUN, NEUT, SETARR)

Create two functions to avoid method ambiguities for a lazy set operation with respect to the neutral element set type and the array set type.

Input

• FUN – function name
• NEUT – set type for neutral element
• SETARR – array set type

Output

A quoted expression containing the function definitions.

Examples

@array_neutral(MinkowskiSum, N, ARR) creates the following functions as quoted expressions:

• MinkowskiSum(N, ARR) = ARR
• MinkowskiSum(ARR, N) = ARR

@array_absorbing(FUN, ABS, SETARR)

Create two functions to avoid method ambiguities for a lazy set operation with respect to the absorbing element set type and the array set type.

Input

• FUN – function name
• ABS – set type for absorbing element
• SETARR – array set type

Output

A quoted expression containing the function definitions.

Examples

@array_absorbing(MinkowskiSum, ABS, ARR) creates the following functions as quoted expressions:

• MinkowskiSum(ABS, ARR) = ABS
• MinkowskiSum(ARR, ABS) = ABS

CachedPair{N}

A mutable pair of an index and a vector.

Fields

• idx – index
• vec – vector

StrictlyIncreasingIndices

Iterator over the vectors of m strictly increasing indices from 1 to n.

Fields

• n – size of the index domain
• m – number of indices to choose (resp. length of the vectors)

Notes

The vectors are modified in-place.

The iterator ranges over $\binom{n}{m}$ (n choose m) possible vectors.

This implementation results in a lexicographic order with the last index growing first.

Examples

julia> for v in LazySets.StrictlyIncreasingIndices(4, 2)
           println(v)
       end
[1, 2]
[1, 3]
[1, 4]
[2, 3]
[2, 4]
[3, 4]

subtypes(interface, concrete::Bool)

Return the concrete subtypes of a given interface.

Input

• interface – an abstract type, usually a set interface
• concrete – if true, seek further the inner abstract subtypes of the given interface, otherwise return only the direct subtypes of interface

Output

A list with the subtypes of the abstract type interface, sorted alphabetically.

Examples

Consider the AbstractPolytope interface. If we include the abstract subtypes of this interface,

julia> using LazySets: subtypes

julia> subtypes(AbstractPolytope, false)
4-element Array{Any,1}:
 AbstractCentrallySymmetricPolytope
 AbstractPolygon
 HPolytope
 VPolytope

We can use this function to obtain the concrete subtypes of AbstractCentrallySymmetricPolytope and AbstractPolygon (further until all concrete types are obtained), using the concrete flag:

julia> subtypes(AbstractPolytope, true)
16-element Array{Type,1}:
 Ball1
 BallInf
 HParallelotope
 HPolygon
 HPolygonOpt
 HPolytope
 Hyperrectangle
 Interval
 LineSegment
 RotatedHyperrectangle
 Singleton
 SymmetricIntervalHull
 VPolygon
 VPolytope
 ZeroSet
 Zonotope

implementing_sets(op::Function;
                  signature::Tuple{Vector{Type}, Int}=(Type[], 1),
                  type_args=Float64, binary::Bool=false)

Compute a dictionary containing information about availability of (unary or binary) concrete set operations.

Input

• op – set operation (respectively its Function object)
• signature – (optional, default: Type[]) the type signature of the function without the LazySet type(s) (see also the index option and the Examples section below)
• index – (optional, default: 1) index of the set type in the signature in the unary case (see the binary option)
• type_args – (optional, default: Float64) type arguments added to the LazySet(s) when searching for available methods; valid inputs are a type or nothing, and in the unary case (see the binary option) it can also be a list of types
• binary – (optional, default: false) flag indicating whether op is a binary function (true) or a unary function (false)

Output

A dictionary with three keys each mapping to a list:

• "available" – This list contains all set types such that there exists an implementation of op.
• "missing" – This list contains all set types such that there does not exist an implementation of op. Note that this is the complement of the "available" list.
• "specific" – This list contains all set types such that there exists a type-specific implementation. Note that those set types also occur in the "available" list.

In the unary case, the lists contain set types. In the binary case, the lists contain pairs of set types.

Examples

julia> using LazySets: implementing_sets

julia> dict = implementing_sets(tovrep);

julia> dict["available"]  # tovrep is only available for polyhedral set types
6-element Array{Type,1}:
 HPolygon
 HPolygonOpt
 HPolyhedron
 HPolytope
 VPolygon
 VPolytope

julia> dict = implementing_sets(σ; signature=Type[AbstractVector{Float64}], index=2);

julia> isempty(dict["missing"])  # every set type implements function σ
true

julia> N = Rational{Int};  # restriction of the number type

julia> dict = implementing_sets(σ; signature=Type[AbstractVector{N}], index=2, type_args=N);

julia> dict["missing"]  # some set types are not available with number type N
3-element Array{Type,1}:
 Ball2
 Ballp
 Ellipsoid

julia> dict = LazySets.implementing_sets(convex_hull; binary=true);  # binary case

julia> (HPolytope, HPolytope) ∈ dict["available"]  # dict contains pairs now
true

read_gen(filename::String)

Read a sequence of polygons stored in vertex representation (gen format).

Input

• filename – path of the file containing the polygons

Output

A list of polygons in vertex representation.

Notes

The x and y coordinates of each vertex should be separated by an empty space, and polygons are separated by empty lines (even the last polygon). For example:

1.01 1.01
0.99 1.01
0.99 0.99
1.01 0.99

0.908463 1.31047
0.873089 1.31047
0.873089 1.28452
0.908463 1.28452

This is parsed as

2-element Array{VPolygon{Float64,Array{Float64,1}},1}:
 VPolygon{Float64,Array{Float64,1}}([[1.01, 1.01], [0.99, 1.01], [0.99, 0.99], [1.01, 0.99]])
 VPolygon{Float64,Array{Float64,1}}([[0.908463, 1.31047], [0.873089, 1.31047], [0.873089, 1.28452], [0.908463, 1.28452]])

_sample_unit_nsphere_muller!(D::Vector{Vector{N}}, n::Int, p::Int;
                             rng::AbstractRNG=GLOBAL_RNG,
                             seed::Union{Int, Nothing}=nothing) where {N}

Draw samples from a uniform distribution on an $n$-dimensional unit sphere using Muller's method.

Input

• D – output, vector of points
• n – dimension of the sphere
• p – number of random samples
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A vector of nsamples vectors.

Algorithm

This function implements Muller's method of normalised Gaussians [1] to uniformly sample over the $n$-dimensional sphere $S^n$ (which is the bounding surface of the $n$-dimensional unit ball).

Given $n$ canonical Gaussian random variables $Z₁, Z₂, …, Z_n$, the distribution of the vectors

\[\dfrac{1}{α}\left(z₁, z₂, …, z_n\right)^T,\]

where $α := \sqrt{z₁² + z₂² + … + z_n²}$, is uniform over $S^n$.

[1] Muller, Mervin E. A note on a method for generating points uniformly on n-dimensional spheres. Communications of the ACM 2.4 (1959): 19-20.

_sample_unit_nball_muller!(D::Vector{Vector{N}}, n::Int, p::Int;
                           rng::AbstractRNG=GLOBAL_RNG,
                           seed::Union{Int, Nothing}=nothing) where {N}

Draw samples from a uniform distribution on an $n$-dimensional unit ball using Muller's method.

Input

• D – output, vector of points
• n – dimension of the ball
• p – number of random samples
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A vector of nsamples vectors.

Algorithm

This function implements Muller's method of normalised Gaussians [1] to sample from the interior of the ball.

Given $n$ Gaussian random variables $Z₁, Z₂, …, Z_n$ and a uniformly distributed random variable $r$ with support in $[0, 1]$, the distribution of the vectors

\[\dfrac{r^{1/n}}{α} \left(z₁, z₂, …, z_n\right)^T,\]

where $α := \sqrt{z₁² + z₂² + … + z_n²}$, is uniform over the $n$-dimensional unit ball.

[1] Muller, Mervin E. A note on a method for generating points uniformly on n-dimensional spheres. Communications of the ACM 2.4 (1959): 19-20.

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

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.

sample(X::LazySet{N}, num_samples::Int;
       [sampler]=_default_sampler(X),
       [rng]::AbstractRNG=GLOBAL_RNG,
       [seed]::Union{Int, Nothing}=nothing,
       [include_vertices]=false,
       [VN]=Vector{N}) where {N}

Sampling of an arbitrary bounded set X.

Input

• X – (bounded) set to be sampled
• num_samples – number of random samples
• sampler – (optional, default: _default_sampler(X)) the sampler used; falls back to RejectionSampler or UniformSampler depending on the type of X
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding
• include_vertices – (optional, default: false) option to include the vertices
• VN – (optional, default: Vector{N}) vector type of the sampled points

Output

A vector of num_samples vectors. If num_samples is not passed, the result is just one sample (not wrapped in a vector).

Algorithm

See the documentation of the respective Sampler.

Notes

If include_vertices == true, we include all vertices computed with vertices. Alternatively if a number $k$ is passed, we plot the first $k$ vertices returned by vertices.

_sample!(D::Vector{VN},
         sampler::RejectionSampler;
         rng::AbstractRNG=GLOBAL_RNG,
         seed::Union{Int, Nothing}=nothing) where {N, VN<:AbstractVector{N}}

Sample points using rejection sampling.

Input

• D – output, vector of points
• sampler – Sampler from which the points are sampled
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A vector of num_samples vectors.

_sample!(D::Vector{VN},
         sampler::UniformSampler{<:LineSegment, <:DefaultUniform};
         rng::AbstractRNG=GLOBAL_RNG,
         seed::Union{Int, Nothing}=nothing) where {N, VN<:AbstractVector{N}}

Sample points of a line segment using uniform sampling in-place.

Input

• D – output, vector of points
• sampler – Sampler from which the points are sampled
• rng – (optional, default: GLOBAL_RNG) random number generator
• seed – (optional, default: nothing) seed for reseeding

Output

A vector of num_samples vectors, where num_samples is the length of D.

Sampler

Abstract type for defining new sample methods.

Notes

All subtypes should implement a _sample! method.

RejectionSampler{S<:LazySet, D} <: Sampler

Type used for rejection sampling of an arbitrary LazySet X.

Fields

• X – (bounded) set to be sampled
• box_approx – Distribution from which the sample is drawn

Algorithm

Draw a sample $x$ from a uniform distribution of a box-overapproximation of the original set $X$ in all $n$ dimensions. The function rejects a drawn sample $x$ and redraws as long as the sample is not contained in the original set $X$, i.e., $x ∉ X$.

UniformSampler{S<:LazySet, D} <: Sampler

Type used for uniform sampling of an arbitrary LazySet X.

Fields

• X – set to be sampled

Algorithm

Draw a sample $x$ from a uniform distribution of the original set $X$.

PolytopeSampler{S<:LazySet, D} <: Sampler

Type used for sampling of a convex polytope X.

Fields

• X – convex polytope to be sampled

Algorithm

Choose a random convex combination of the vertices of X.

Let $V = \{v_i\}_i$ denote the set of vertices of X. Then any point $p \in \mathbb{R}^n$ of the convex polytope $X$ is a convex combination of its vertices, i.e., $p = \sum_{i} v_i α_i$ for some (non-negative) coefficients $\{α_i\}_i$ that add up to 1. The algorithm chooses a random convex combination (the $α_i$). To produce such combination we apply the finite difference operator on a sorted uniform sample over $[0, 1]$; the method can be found in [1] and [2].

Notes

The sampling is not uniform - points in the center of the polytope are more likely to be sampled.

References

[1] Rubin, Donald B. The bayesian bootstrap. The annals of statistics (1981): 130-134.

[2] https://cs.stackexchange.com/questions/3227/uniform-sampling-from-a-simplex/3229

volume(H::AbstractHyperrectangle)

Return the volume of a hyperrectangular set.

Input

• H – hyperrectangular set

Output

The volume of $H$.

Algorithm

The volume of the $n$-dimensional hyperrectangle $H$ with vector radius $r$ is $2ⁿ ∏ᵢ rᵢ$ where $rᵢ$ denotes the $i$-th component of $r$.

volume(B::Ball2)

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

Input

• B – ball in the 2-norm

Output

The volume of $B$.

Algorithm

This function implements the well-known formula for the volume of an n-dimensional ball using factorials. For details see the wikipedia article Volume of an n-ball.

_vec(vars::NTuple{L, Union{<:Num, <:Vector{Num}}}) where {L}

Transform a tuple of operations into one vector of operations.

Input

• vars – tuple where each element is either variable-like (Num) or a vector of variables (Vector{Num})

Output

A vector of Operation obtained by concatenating each tuple component.

Examples

julia> vars = @variables x[1:2] y
(Num[x₁, x₂], y)

julia> LazySets._vec(vars)
3-element Array{Num,1}:
 x₁
 x₂
 y