From e7202ec5e53773d94cb99c48992a8d908c1f0c92 Mon Sep 17 00:00:00 2001 From: schillic Date: Fri, 12 Jul 2024 18:03:27 +0200 Subject: [PATCH 1/6] move HParallelotope.jl to subfolder --- src/LazySets.jl | 3 ++- .../HParallelotopeModule.jl} | 0 2 files changed, 2 insertions(+), 1 deletion(-) rename src/Sets/{HParallelotope.jl => HParallelotope/HParallelotopeModule.jl} (100%) diff --git a/src/LazySets.jl b/src/LazySets.jl index 8d8161dcd9..b3c0487851 100644 --- a/src/LazySets.jl +++ b/src/LazySets.jl @@ -119,7 +119,8 @@ include("Sets/Ellipsoid/EllipsoidModule.jl") include("Sets/DensePolynomialZonotope/DensePolynomialZonotopeModule.jl") @reexport using ..DensePolynomialZonotopeModule: DensePolynomialZonotope -include("Sets/HParallelotope.jl") +include("Sets/HParallelotope/HParallelotopeModule.jl") + include("Sets/HPolygon.jl") include("Sets/HPolygonOpt.jl") include("Sets/HPolytope.jl") diff --git a/src/Sets/HParallelotope.jl b/src/Sets/HParallelotope/HParallelotopeModule.jl similarity index 100% rename from src/Sets/HParallelotope.jl rename to src/Sets/HParallelotope/HParallelotopeModule.jl From 236370561cc7cf4c74ef74a0d9010fe66c45e885 Mon Sep 17 00:00:00 2001 From: schillic Date: Fri, 12 Jul 2024 19:24:42 +0200 Subject: [PATCH 2/6] add HParallelotopeModule --- docs/src/lib/sets/HParallelotope.md | 43 +++++++++++-------- src/LazySets.jl | 5 +++ .../HParallelotope/HParallelotopeModule.jl | 21 +++++++-- 3 files changed, 49 insertions(+), 20 deletions(-) diff --git a/docs/src/lib/sets/HParallelotope.md b/docs/src/lib/sets/HParallelotope.md index fbe32ce94a..db33966f45 100644 --- a/docs/src/lib/sets/HParallelotope.md +++ b/docs/src/lib/sets/HParallelotope.md @@ -1,49 +1,58 @@ ```@meta -CurrentModule = LazySets +CurrentModule = LazySets.HParallelotopeModule ``` # [HParallelotope](@id def_HParallelotope) ```@docs HParallelotope -directions(::HParallelotope) -offset(::HParallelotope) -dim(::HParallelotope) +``` + +## Operations + +```@docs base_vertex(::HParallelotope) -extremal_vertices(::HParallelotope{N, VN}) where {N, VN} center(::HParallelotope) -genmat(::HParallelotope) -generators(::HParallelotope) constraints_list(::HParallelotope) +dim(::HParallelotope) +directions(::HParallelotope) +extremal_vertices(::HParallelotope{N, VN}) where {N, VN} +generators(::HParallelotope) +genmat(::HParallelotope) +offset(::HParallelotope) rand(::Type{HParallelotope}) volume(::HParallelotope) ``` +```@meta +CurrentModule = LazySets +``` + Inherited from [`LazySet`](@ref): -* [`low`](@ref low(::LazySet)) +* [`diameter`](@ref diameter(::LazySet, ::Real)) * [`high`](@ref high(::LazySet)) +* [`low`](@ref low(::LazySet)) * [`norm`](@ref norm(::LazySet, ::Real)) * [`radius`](@ref radius(::LazySet, ::Real)) -* [`diameter`](@ref diameter(::LazySet, ::Real)) -* [`singleton_list`](@ref singleton_list(::LazySet)) * [`rectify`](@ref rectify(::LazySet)) +* [`singleton_list`](@ref singleton_list(::LazySet)) Inherited from [`AbstractPolytope`](@ref): * [`isbounded`](@ref isbounded(::AbstractPolytope)) Inherited from [`AbstractCentrallySymmetricPolytope`](@ref): -* [`isuniversal`](@ref isuniversal(::AbstractCentrallySymmetricPolytope{N}, ::Bool=false) where {N}) * [`an_element`](@ref an_element(::AbstractCentrallySymmetricPolytope)) * [`extrema`](@ref extrema(::AbstractCentrallySymmetricPolytope)) * [`extrema`](@ref extrema(::AbstractCentrallySymmetricPolytope, ::Int)) +* [`isuniversal`](@ref isuniversal(::AbstractCentrallySymmetricPolytope{N}, ::Bool=false) where {N}) Inherited from [`AbstractZonotope`](@ref): -* [`ρ`](@ref ρ(::AbstractVector, ::AbstractZonotope)) -* [`σ`](@ref σ(::AbstractVector, ::AbstractZonotope)) -* [`∈`](@ref ∈(::AbstractVector, ::AbstractZonotope)) -* [`linear_map`](@ref linear_map(::AbstractMatrix, ::AbstractZonotope)) -* [`translate`](@ref translate(::AbstractZonotope, ::AbstractVector)) -* [`vertices_list`](@ref vertices_list(::AbstractZonotope)) * [`isempty`](@ref isempty(::AbstractZonotope)) * [`order`](@ref order(::AbstractZonotope)) * [`reflect`](@ref reflect(::AbstractZonotope)) +* [`vertices_list`](@ref vertices_list(::AbstractZonotope)) +* [`linear_map`](@ref linear_map(::AbstractMatrix, ::AbstractZonotope)) +* [`∈`](@ref ∈(::AbstractVector, ::AbstractZonotope)) +* [`ρ`](@ref ρ(::AbstractVector, ::AbstractZonotope)) +* [`σ`](@ref σ(::AbstractVector, ::AbstractZonotope)) +* [`translate`](@ref translate(::AbstractZonotope, ::AbstractVector)) diff --git a/src/LazySets.jl b/src/LazySets.jl index b3c0487851..165abc2712 100644 --- a/src/LazySets.jl +++ b/src/LazySets.jl @@ -120,6 +120,11 @@ include("Sets/DensePolynomialZonotope/DensePolynomialZonotopeModule.jl") @reexport using ..DensePolynomialZonotopeModule: DensePolynomialZonotope include("Sets/HParallelotope/HParallelotopeModule.jl") +@reexport using ..HParallelotopeModule: HParallelotope, + directions, + base_vertex, + extremal_vertices, + offset include("Sets/HPolygon.jl") include("Sets/HPolygonOpt.jl") diff --git a/src/Sets/HParallelotope/HParallelotopeModule.jl b/src/Sets/HParallelotope/HParallelotopeModule.jl index 1a3ef4b169..c886f454c8 100644 --- a/src/Sets/HParallelotope/HParallelotopeModule.jl +++ b/src/Sets/HParallelotope/HParallelotopeModule.jl @@ -1,10 +1,23 @@ +module HParallelotopeModule + +using Reexport + +using ..LazySets: AbstractZonotope, HalfSpace, HPolyhedron, generators_fallback +using Random: AbstractRNG, GLOBAL_RNG +using ReachabilityBase.Arrays: to_negative_vector +using ReachabilityBase.Distribution: reseed! +using LinearAlgebra: checksquare, det + +@reexport import ..API: center, constraints_list, dim, isoperationtype, rand, + volume +@reexport import ..LazySets: generators, genmat +@reexport using ..API + export HParallelotope, directions, - offset, base_vertex, extremal_vertices, - genmat, - generators + offset """ HParallelotope{N, VN<:AbstractVector{N}, MN<:AbstractMatrix{N}} <: AbstractZonotope{N} @@ -423,3 +436,5 @@ function volume(P::HParallelotope) n = dim(P) return 2^n * abs(det(genmat(P))) end + +end # module From 6ea488071d69d2033d8496bfddbb495b82ccb423 Mon Sep 17 00:00:00 2001 From: schillic Date: Fri, 12 Jul 2024 19:32:04 +0200 Subject: [PATCH 3/6] outsource struct and operations to files --- src/Sets/HParallelotope/HParallelotope.jl | 78 ++++ .../HParallelotope/HParallelotopeModule.jl | 431 +----------------- src/Sets/HParallelotope/base_vertex.jl | 29 ++ src/Sets/HParallelotope/center.jl | 31 ++ src/Sets/HParallelotope/constraints_list.jl | 36 ++ src/Sets/HParallelotope/dim.jl | 16 + src/Sets/HParallelotope/directions.jl | 18 + src/Sets/HParallelotope/extremal_vertices.jl | 45 ++ src/Sets/HParallelotope/generators.jl | 17 + src/Sets/HParallelotope/genmat.jl | 30 ++ src/Sets/HParallelotope/isoperationtype.jl | 1 + src/Sets/HParallelotope/offset.jl | 17 + src/Sets/HParallelotope/rand.jl | 57 +++ src/Sets/HParallelotope/volume.jl | 26 ++ 14 files changed, 415 insertions(+), 417 deletions(-) create mode 100644 src/Sets/HParallelotope/HParallelotope.jl create mode 100644 src/Sets/HParallelotope/base_vertex.jl create mode 100644 src/Sets/HParallelotope/center.jl create mode 100644 src/Sets/HParallelotope/constraints_list.jl create mode 100644 src/Sets/HParallelotope/dim.jl create mode 100644 src/Sets/HParallelotope/directions.jl create mode 100644 src/Sets/HParallelotope/extremal_vertices.jl create mode 100644 src/Sets/HParallelotope/generators.jl create mode 100644 src/Sets/HParallelotope/genmat.jl create mode 100644 src/Sets/HParallelotope/isoperationtype.jl create mode 100644 src/Sets/HParallelotope/offset.jl create mode 100644 src/Sets/HParallelotope/rand.jl create mode 100644 src/Sets/HParallelotope/volume.jl diff --git a/src/Sets/HParallelotope/HParallelotope.jl b/src/Sets/HParallelotope/HParallelotope.jl new file mode 100644 index 0000000000..55aba9239c --- /dev/null +++ b/src/Sets/HParallelotope/HParallelotope.jl @@ -0,0 +1,78 @@ +""" + HParallelotope{N, VN<:AbstractVector{N}, MN<:AbstractMatrix{N}} <: AbstractZonotope{N} + +Type that represents a parallelotope in constraint form. + +### Fields + +- `directions` -- square matrix where each row is the direction of two parallel + constraints +- `offset` -- vector where each element is the offset of the corresponding + constraint + +### Notes + +Parallelotopes are centrally symmetric convex polytopes in ``ℝ^n`` +having ``2n`` pairwise parallel constraints. Every parallelotope is a zonotope. +As such, parallelotopes can be represented in constraint form or in generator +form. The `HParallelotope` type represents parallelotopes in constraint form. + +Let ``D ∈ ℝ^{n × n}`` be a matrix and let ``c ∈ ℝ^{2n}`` be +a vector. The parallelotope ``P ⊂ ℝ^n`` generated by the directions +matrix ``D`` and the offset vector ``c`` is given by the set of points +``x ∈ ℝ^`` such that: + +```math + D_i ⋅ x ≤ c_{i},\\text{ and } -D_i ⋅ x ≤ c_{n+i} +``` +for ``i = 1, …, n``. Here ``D_i`` represents the ``i``-th row of ``D`` and +``c_i`` the ``i``-th component of ``c``. + +Note that, although representing a zonotopic set, an `HParallelotope` can be +empty or unbounded if the constraints are unsuitably chosen. This may cause +problems with default methods because the library assumes that zonotopic sets +are non-empty and bounded. Thus such instances are considered illegal. The +default constructor thus checks these conditions, which can be deactivated by +passing the argument `check_consistency=false`. + +For details as well as applications of parallelotopes in reachability analysis +we refer to [1] and [2]. For conversions between set representations we refer to +[3]. + +### References + +[1] Tommaso Dreossi, Thao Dang, and Carla Piazza. *Reachability computation for +polynomial dynamical systems.* Formal Methods in System Design 50.1 (2017): 1-38. + +[2] Tommaso Dreossi, Thao Dang, and Carla Piazza. *Parallelotope bundles for +polynomial reachability.* Proceedings of the 19th International Conference on +Hybrid Systems: Computation and Control. ACM, 2016. + +[3] Matthias Althoff, Olaf Stursberg, and Martin Buss. *Computing reachable sets +of hybrid systems using a combination of zonotopes and polytopes.* Nonlinear +analysis: hybrid systems 4.2 (2010): 233-249. +""" +struct HParallelotope{N,VN<:AbstractVector{N},MN<:AbstractMatrix{N}} <: AbstractZonotope{N} + directions::MN + offset::VN + + # default constructor with dimension and consistency check + function HParallelotope(D::MN, c::VN; + check_consistency::Bool=true) where {N,VN<:AbstractVector{N}, + MN<:AbstractMatrix{N}} + @assert length(c) == 2 * checksquare(D) "the length of the offset " * + "vector should be twice the size of the directions matrix, " * + "but they have sizes $(length(c)) and $(size(D)) respectively" + + if check_consistency + P = HPolyhedron(_constraints_list_hparallelotope(D, c, N, VN)) + if isempty(P) + throw(ArgumentError("the constraints are contradictory")) + elseif !isbounded(P) + throw(ArgumentError("the constraints are not bounding")) + end + end + + return new{N,VN,MN}(D, c) + end +end diff --git a/src/Sets/HParallelotope/HParallelotopeModule.jl b/src/Sets/HParallelotope/HParallelotopeModule.jl index c886f454c8..aeb91e8a97 100644 --- a/src/Sets/HParallelotope/HParallelotopeModule.jl +++ b/src/Sets/HParallelotope/HParallelotopeModule.jl @@ -19,422 +19,19 @@ export HParallelotope, extremal_vertices, offset -""" - HParallelotope{N, VN<:AbstractVector{N}, MN<:AbstractMatrix{N}} <: AbstractZonotope{N} - -Type that represents a parallelotope in constraint form. - -### Fields - -- `directions` -- square matrix where each row is the direction of two parallel - constraints -- `offset` -- vector where each element is the offset of the corresponding - constraint - -### Notes - -Parallelotopes are centrally symmetric convex polytopes in ``ℝ^n`` -having ``2n`` pairwise parallel constraints. Every parallelotope is a zonotope. -As such, parallelotopes can be represented in constraint form or in generator -form. The `HParallelotope` type represents parallelotopes in constraint form. - -Let ``D ∈ ℝ^{n × n}`` be a matrix and let ``c ∈ ℝ^{2n}`` be -a vector. The parallelotope ``P ⊂ ℝ^n`` generated by the directions -matrix ``D`` and the offset vector ``c`` is given by the set of points -``x ∈ ℝ^`` such that: - -```math - D_i ⋅ x ≤ c_{i},\\text{ and } -D_i ⋅ x ≤ c_{n+i} -``` -for ``i = 1, …, n``. Here ``D_i`` represents the ``i``-th row of ``D`` and -``c_i`` the ``i``-th component of ``c``. - -Note that, although representing a zonotopic set, an `HParallelotope` can be -empty or unbounded if the constraints are unsuitably chosen. This may cause -problems with default methods because the library assumes that zonotopic sets -are non-empty and bounded. Thus such instances are considered illegal. The -default constructor thus checks these conditions, which can be deactivated by -passing the argument `check_consistency=false`. - -For details as well as applications of parallelotopes in reachability analysis -we refer to [1] and [2]. For conversions between set representations we refer to -[3]. - -### References - -[1] Tommaso Dreossi, Thao Dang, and Carla Piazza. *Reachability computation for -polynomial dynamical systems.* Formal Methods in System Design 50.1 (2017): 1-38. - -[2] Tommaso Dreossi, Thao Dang, and Carla Piazza. *Parallelotope bundles for -polynomial reachability.* Proceedings of the 19th International Conference on -Hybrid Systems: Computation and Control. ACM, 2016. - -[3] Matthias Althoff, Olaf Stursberg, and Martin Buss. *Computing reachable sets -of hybrid systems using a combination of zonotopes and polytopes.* Nonlinear -analysis: hybrid systems 4.2 (2010): 233-249. -""" -struct HParallelotope{N,VN<:AbstractVector{N},MN<:AbstractMatrix{N}} <: AbstractZonotope{N} - directions::MN - offset::VN - - # default constructor with dimension and consistency check - function HParallelotope(D::MN, c::VN; - check_consistency::Bool=true) where {N,VN<:AbstractVector{N}, - MN<:AbstractMatrix{N}} - @assert length(c) == 2 * checksquare(D) "the length of the offset " * - "vector should be twice the size of the directions matrix, " * - "but they have sizes $(length(c)) and $(size(D)) respectively" - - if check_consistency - P = HPolyhedron(_constraints_list_hparallelotope(D, c, N, VN)) - if isempty(P) - throw(ArgumentError("the constraints are contradictory")) - elseif !isbounded(P) - throw(ArgumentError("the constraints are not bounding")) - end - end - - return new{N,VN,MN}(D, c) - end -end - -isoperationtype(::Type{<:HParallelotope}) = false - -# ================= -# Getter functions -# ================= - -""" - directions(P::HParallelotope) - -Return the directions matrix of a parallelotope in constraint representation. - -### Input - -- `P` -- parallelotope in constraint representation - -### Output - -A matrix where each row represents a direction of the parallelotope. -The negated directions `-D_i` are implicit (see [`HParallelotope`](@ref) for -details). -""" -function directions(P::HParallelotope) - return P.directions -end - -""" - offset(P::HParallelotope) - -Return the offsets of a parallelotope in constraint representation. - -### Input - -- `P` -- parallelotope in constraint representation - -### Output - -A vector with the ``2n`` offsets of the parallelotope, where ``n`` is the -dimension of `P`. -""" -function offset(P::HParallelotope) - return P.offset -end - -""" - dim(P::HParallelotope) - -Return the dimension of a parallelotope in constraint representation. - -### Input - -- `P` -- parallelotope in constraint representation - -### Output - -The ambient dimension of the parallelotope. -""" -function dim(P::HParallelotope) - return size(P.directions, 1) -end - -""" - base_vertex(P::HParallelotope) - -Compute the base vertex of a parallelotope in constraint representation. - -### Input - -- `P` -- parallelotope in constraint representation - -### Output - -The base vertex of `P`. - -### Algorithm - -Intuitively, the base vertex is the point from which we get the relative -positions of all the other points. -The base vertex can be computed as the solution of the ``n``-dimensional linear -system ``D_i x = c_{n+i}`` for ``i = 1, …, n``, see [1, Section 3.2.1]. - -[1] Dreossi, Tommaso, Thao Dang, and Carla Piazza. *Reachability computation for -polynomial dynamical systems.* Formal Methods in System Design 50.1 (2017): 1-38. -""" -function base_vertex(P::HParallelotope) - D, c = P.directions, P.offset - n = dim(P) - v = to_negative_vector(view(c, (n + 1):(2n))) # converts to a dense vector as well - return D \ v -end - -""" - extremal_vertices(P::HParallelotope{N, VN}) where {N, VN} - -Compute the extremal vertices with respect to the base vertex of a parallelotope -in constraint representation. - -### Input - -- `P` -- parallelotope in constraint representation - -### Output - -The list of vertices connected to the base vertex of ``P``. - -### Notes - -Let ``P`` be a parallelotope in constraint representation with directions matrix -``D`` and offset vector ``c``. The *extremal vertices* of ``P`` with respect to -its base vertex ``q`` are those vertices of ``P`` that have an edge in common -with ``q``. - -### Algorithm - -The extremal vertices can be computed as the solution of the ``n``-dimensional -linear systems of equations ``D x = v_i`` where for each ``i = 1, …, n``, -``v_i = [-c_{n+1}, …, c_i, …, -c_{2n}]``. - -We refer to [1, Section 3.2.1] for details. - -[1] Tommaso Dreossi, Thao Dang, and Carla Piazza. *Reachability computation for -polynomial dynamical systems.* Formal Methods in System Design 50.1 (2017): 1-38. -""" -function extremal_vertices(P::HParallelotope{N,VN}) where {N,VN} - D, c = P.directions, P.offset - n = dim(P) - v = to_negative_vector(view(c, (n + 1):(2n))) - vertices = Vector{VN}(undef, n) - h = copy(v) - @inbounds for i in 1:n - h[i] = c[i] - vertices[i] = D \ h - h[i] = v[i] - end - return vertices -end - -""" - center(P::HParallelotope) - -Return the center of a parallelotope in constraint representation. - -### Input - -- `P` -- parallelotope in constraint representation - -### Output - -The center of the parallelotope. - -### Algorithm - -Let ``P`` be a parallelotope with base vertex ``q`` and list of extremal -vertices with respect to ``q`` given by the set ``\\{v_i\\}`` for -``i = 1, …, n``. Then the center is located at - -```math - c = q + ∑_{i=1}^n \\frac{v_i - q}{2} = q (1 - \\frac{2}) + \\frac{s}{2}, -``` -where ``s := ∑_{i=1}^n v_i`` is the sum of extremal vertices. -""" -function center(P::HParallelotope) - n = dim(P) - q = base_vertex(P) - E = extremal_vertices(P) - s = sum(E) - return q * (1 - n / 2) + s / 2 -end - -""" - genmat(P::HParallelotope) - -Return the generator matrix of a parallelotope in constraint representation. - -### Input - -- `P` -- parallelotope in constraint representation - -### Output - -A matrix where each column represents one generator of the parallelotope `P`. - -### Algorithm - -Let ``P`` be a parallelotope with base vertex ``q`` and list of extremal -vertices with respect to ``q`` given by the set ``\\{v_i\\}`` for -``i = 1, …, n``. Then, the ``i``-th generator of ``P``, represented as the -``i``-th column vector ``G[:, i]``, is given by: - -```math - G[:, i] = \\frac{v_i - q}{2} -``` -for ``i = 1, …, n``. -""" -function genmat(P::HParallelotope) - q = base_vertex(P) - E = extremal_vertices(P) - return 1 / 2 * reduce(hcat, E) .- q / 2 -end - -""" - generators(P::HParallelotope) - -Return an iterator over the generators of a parallelotope in constraint -representation. - -### Input - -- `P` -- parallelotope in constraint representation - -### Output - -An iterator over the generators of `P`. -""" -function generators(P::HParallelotope) - return generators_fallback(P) -end - -""" - constraints_list(P::HParallelotope) - -Return the list of constraints of a parallelotope in constraint representation. - -### Input - -- `P` -- parallelotope in constraint representation - -### Output - -The list of constraints of `P`. -""" -function constraints_list(P::HParallelotope) - D, c = P.directions, P.offset - N, VN = _parameters(P) - return _constraints_list_hparallelotope(D, c, N, VN) -end - -function _constraints_list_hparallelotope(D, c, N, VN) - if isempty(D) - return Vector{HalfSpace{N,VN}}(undef, 0) - end - n = size(D, 1) - clist = Vector{HalfSpace{N,VN}}(undef, 2n) - @inbounds for i in 1:n - clist[i] = HalfSpace(D[i, :], c[i]) - clist[i + n] = HalfSpace(-D[i, :], c[i + n]) - end - return clist -end - -# reason: `Documenter` cannot annotate `constraints_list` with type parameters -function _parameters(P::HParallelotope{N,VN}) where {N,VN} - return (N, VN) -end - -""" - rand(::Type{HParallelotope}; [N]::Type{<:Real}=Float64, [dim]::Int=2, - [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing) - -Create a random parallelotope in constraint representation. - -### Input - -- `HParallelotope` -- type for dispatch -- `N` -- (optional, default: `Float64`) numeric type -- `dim` -- (optional, default: 2) dimension -- `rng` -- (optional, default: `GLOBAL_RNG`) random number generator -- `seed` -- (optional, default: `nothing`) seed for reseeding - -### Output - -A random parallelotope. - -### Notes - -All numbers are normally distributed with mean 0 and standard deviation 1. - -### Algorithm - -The directions matrix and offset vector are created randomly. On average there -is a good chance that this resulting set is empty. We then modify the offset to -ensure non-emptiness. - -There is a chance that the resulting set represents an unbounded set. This -implementation checks for that case and then samples a new set. -""" -function rand(::Type{HParallelotope}; - N::Type{<:Real}=Float64, - dim::Int=2, - rng::AbstractRNG=GLOBAL_RNG, - seed::Union{Int,Nothing}=nothing) - rng = reseed!(rng, seed) - - while true - D = randn(rng, N, dim, dim) - offset = randn(rng, N, 2 * dim) - - # make sure that the set is not empty: - # `offset[i] >= -offset[i-dim]` for all `i ∈ dim+1:2*dim` - @inbounds for i in (dim + 1):(2 * dim) - offset[i] = -offset[i - dim] + abs(offset[i]) - end - - # convert to polyhedron to check boundedness - clist = _constraints_list_hparallelotope(D, offset, N, typeof(offset)) - Q = HPolyhedron(clist) - if isbounded(Q) - return P = HParallelotope(D, offset; check_consistency=false) - end - # set is unbounded; sample a new set in the next iteration - end -end - -""" - volume(P::HParallelotope) - -Return the volume of a parallelotope in constraint representation. - -### Input - -- `P` -- parallelotope in constraint representation - -### Output - -The volume. - -### Algorithm - -The volume of an ``n``-dimensional parallelotope `P` is ``2^n · |\\det(G)|``, -where ``G`` is the generator matrix of `P`. This can be seen as follows: -The generator matrix transforms the ``n``-dimensional hypercube ``[0, 1]^n`` -to a parallelotope of volume ``|\\det(G)|``. Since the representation of a -parallelotope instead transforms the hypercube ``[-1, 1]^n``, this result has to -be doubled for each dimension. -""" -function volume(P::HParallelotope) - n = dim(P) - return 2^n * abs(det(genmat(P))) -end +include("HParallelotope.jl") + +include("base_vertex.jl") +include("center.jl") +include("constraints_list.jl") +include("dim.jl") +include("directions.jl") +include("extremal_vertices.jl") +include("generators.jl") +include("genmat.jl") +include("isoperationtype.jl") +include("offset.jl") +include("rand.jl") +include("volume.jl") end # module diff --git a/src/Sets/HParallelotope/base_vertex.jl b/src/Sets/HParallelotope/base_vertex.jl new file mode 100644 index 0000000000..52d7ae4aa9 --- /dev/null +++ b/src/Sets/HParallelotope/base_vertex.jl @@ -0,0 +1,29 @@ +""" + base_vertex(P::HParallelotope) + +Compute the base vertex of a parallelotope in constraint representation. + +### Input + +- `P` -- parallelotope in constraint representation + +### Output + +The base vertex of `P`. + +### Algorithm + +Intuitively, the base vertex is the point from which we get the relative +positions of all the other points. +The base vertex can be computed as the solution of the ``n``-dimensional linear +system ``D_i x = c_{n+i}`` for ``i = 1, …, n``, see [1, Section 3.2.1]. + +[1] Dreossi, Tommaso, Thao Dang, and Carla Piazza. *Reachability computation for +polynomial dynamical systems.* Formal Methods in System Design 50.1 (2017): 1-38. +""" +function base_vertex(P::HParallelotope) + D, c = P.directions, P.offset + n = dim(P) + v = to_negative_vector(view(c, (n + 1):(2n))) # converts to a dense vector as well + return D \ v +end diff --git a/src/Sets/HParallelotope/center.jl b/src/Sets/HParallelotope/center.jl new file mode 100644 index 0000000000..84422cef68 --- /dev/null +++ b/src/Sets/HParallelotope/center.jl @@ -0,0 +1,31 @@ +""" + center(P::HParallelotope) + +Return the center of a parallelotope in constraint representation. + +### Input + +- `P` -- parallelotope in constraint representation + +### Output + +The center of the parallelotope. + +### Algorithm + +Let ``P`` be a parallelotope with base vertex ``q`` and list of extremal +vertices with respect to ``q`` given by the set ``\\{v_i\\}`` for +``i = 1, …, n``. Then the center is located at + +```math + c = q + ∑_{i=1}^n \\frac{v_i - q}{2} = q (1 - \\frac{2}) + \\frac{s}{2}, +``` +where ``s := ∑_{i=1}^n v_i`` is the sum of extremal vertices. +""" +function center(P::HParallelotope) + n = dim(P) + q = base_vertex(P) + E = extremal_vertices(P) + s = sum(E) + return q * (1 - n / 2) + s / 2 +end diff --git a/src/Sets/HParallelotope/constraints_list.jl b/src/Sets/HParallelotope/constraints_list.jl new file mode 100644 index 0000000000..b9492be33c --- /dev/null +++ b/src/Sets/HParallelotope/constraints_list.jl @@ -0,0 +1,36 @@ +""" + constraints_list(P::HParallelotope) + +Return the list of constraints of a parallelotope in constraint representation. + +### Input + +- `P` -- parallelotope in constraint representation + +### Output + +The list of constraints of `P`. +""" +function constraints_list(P::HParallelotope) + D, c = P.directions, P.offset + N, VN = _parameters(P) + return _constraints_list_hparallelotope(D, c, N, VN) +end + +# reason: `Documenter` cannot annotate `constraints_list` with type parameters +function _parameters(P::HParallelotope{N,VN}) where {N,VN} + return (N, VN) +end + +function _constraints_list_hparallelotope(D, c, N, VN) + if isempty(D) + return Vector{HalfSpace{N,VN}}(undef, 0) + end + n = size(D, 1) + clist = Vector{HalfSpace{N,VN}}(undef, 2n) + @inbounds for i in 1:n + clist[i] = HalfSpace(D[i, :], c[i]) + clist[i + n] = HalfSpace(-D[i, :], c[i + n]) + end + return clist +end diff --git a/src/Sets/HParallelotope/dim.jl b/src/Sets/HParallelotope/dim.jl new file mode 100644 index 0000000000..71a9b83770 --- /dev/null +++ b/src/Sets/HParallelotope/dim.jl @@ -0,0 +1,16 @@ +""" + dim(P::HParallelotope) + +Return the dimension of a parallelotope in constraint representation. + +### Input + +- `P` -- parallelotope in constraint representation + +### Output + +The ambient dimension of the parallelotope. +""" +function dim(P::HParallelotope) + return size(P.directions, 1) +end diff --git a/src/Sets/HParallelotope/directions.jl b/src/Sets/HParallelotope/directions.jl new file mode 100644 index 0000000000..bb2bfa22a6 --- /dev/null +++ b/src/Sets/HParallelotope/directions.jl @@ -0,0 +1,18 @@ +""" + directions(P::HParallelotope) + +Return the directions matrix of a parallelotope in constraint representation. + +### Input + +- `P` -- parallelotope in constraint representation + +### Output + +A matrix where each row represents a direction of the parallelotope. +The negated directions `-D_i` are implicit (see [`HParallelotope`](@ref) for +details). +""" +function directions(P::HParallelotope) + return P.directions +end diff --git a/src/Sets/HParallelotope/extremal_vertices.jl b/src/Sets/HParallelotope/extremal_vertices.jl new file mode 100644 index 0000000000..b975517027 --- /dev/null +++ b/src/Sets/HParallelotope/extremal_vertices.jl @@ -0,0 +1,45 @@ +""" + extremal_vertices(P::HParallelotope{N, VN}) where {N, VN} + +Compute the extremal vertices with respect to the base vertex of a parallelotope +in constraint representation. + +### Input + +- `P` -- parallelotope in constraint representation + +### Output + +The list of vertices connected to the base vertex of ``P``. + +### Notes + +Let ``P`` be a parallelotope in constraint representation with directions matrix +``D`` and offset vector ``c``. The *extremal vertices* of ``P`` with respect to +its base vertex ``q`` are those vertices of ``P`` that have an edge in common +with ``q``. + +### Algorithm + +The extremal vertices can be computed as the solution of the ``n``-dimensional +linear systems of equations ``D x = v_i`` where for each ``i = 1, …, n``, +``v_i = [-c_{n+1}, …, c_i, …, -c_{2n}]``. + +We refer to [1, Section 3.2.1] for details. + +[1] Tommaso Dreossi, Thao Dang, and Carla Piazza. *Reachability computation for +polynomial dynamical systems.* Formal Methods in System Design 50.1 (2017): 1-38. +""" +function extremal_vertices(P::HParallelotope{N,VN}) where {N,VN} + D, c = P.directions, P.offset + n = dim(P) + v = to_negative_vector(view(c, (n + 1):(2n))) + vertices = Vector{VN}(undef, n) + h = copy(v) + @inbounds for i in 1:n + h[i] = c[i] + vertices[i] = D \ h + h[i] = v[i] + end + return vertices +end diff --git a/src/Sets/HParallelotope/generators.jl b/src/Sets/HParallelotope/generators.jl new file mode 100644 index 0000000000..300452744e --- /dev/null +++ b/src/Sets/HParallelotope/generators.jl @@ -0,0 +1,17 @@ +""" + generators(P::HParallelotope) + +Return an iterator over the generators of a parallelotope in constraint +representation. + +### Input + +- `P` -- parallelotope in constraint representation + +### Output + +An iterator over the generators of `P`. +""" +function generators(P::HParallelotope) + return generators_fallback(P) +end diff --git a/src/Sets/HParallelotope/genmat.jl b/src/Sets/HParallelotope/genmat.jl new file mode 100644 index 0000000000..53c6d9d71c --- /dev/null +++ b/src/Sets/HParallelotope/genmat.jl @@ -0,0 +1,30 @@ +""" + genmat(P::HParallelotope) + +Return the generator matrix of a parallelotope in constraint representation. + +### Input + +- `P` -- parallelotope in constraint representation + +### Output + +A matrix where each column represents one generator of the parallelotope `P`. + +### Algorithm + +Let ``P`` be a parallelotope with base vertex ``q`` and list of extremal +vertices with respect to ``q`` given by the set ``\\{v_i\\}`` for +``i = 1, …, n``. Then, the ``i``-th generator of ``P``, represented as the +``i``-th column vector ``G[:, i]``, is given by: + +```math + G[:, i] = \\frac{v_i - q}{2} +``` +for ``i = 1, …, n``. +""" +function genmat(P::HParallelotope) + q = base_vertex(P) + E = extremal_vertices(P) + return 1 / 2 * reduce(hcat, E) .- q / 2 +end diff --git a/src/Sets/HParallelotope/isoperationtype.jl b/src/Sets/HParallelotope/isoperationtype.jl new file mode 100644 index 0000000000..43571af362 --- /dev/null +++ b/src/Sets/HParallelotope/isoperationtype.jl @@ -0,0 +1 @@ +isoperationtype(::Type{<:HParallelotope}) = false diff --git a/src/Sets/HParallelotope/offset.jl b/src/Sets/HParallelotope/offset.jl new file mode 100644 index 0000000000..b8c1c1fb80 --- /dev/null +++ b/src/Sets/HParallelotope/offset.jl @@ -0,0 +1,17 @@ +""" + offset(P::HParallelotope) + +Return the offsets of a parallelotope in constraint representation. + +### Input + +- `P` -- parallelotope in constraint representation + +### Output + +A vector with the ``2n`` offsets of the parallelotope, where ``n`` is the +dimension of `P`. +""" +function offset(P::HParallelotope) + return P.offset +end diff --git a/src/Sets/HParallelotope/rand.jl b/src/Sets/HParallelotope/rand.jl new file mode 100644 index 0000000000..0a0a6b673f --- /dev/null +++ b/src/Sets/HParallelotope/rand.jl @@ -0,0 +1,57 @@ +""" + rand(::Type{HParallelotope}; [N]::Type{<:Real}=Float64, [dim]::Int=2, + [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing) + +Create a random parallelotope in constraint representation. + +### Input + +- `HParallelotope` -- type for dispatch +- `N` -- (optional, default: `Float64`) numeric type +- `dim` -- (optional, default: 2) dimension +- `rng` -- (optional, default: `GLOBAL_RNG`) random number generator +- `seed` -- (optional, default: `nothing`) seed for reseeding + +### Output + +A random parallelotope. + +### Notes + +All numbers are normally distributed with mean 0 and standard deviation 1. + +### Algorithm + +The directions matrix and offset vector are created randomly. On average there +is a good chance that this resulting set is empty. We then modify the offset to +ensure non-emptiness. + +There is a chance that the resulting set represents an unbounded set. This +implementation checks for that case and then samples a new set. +""" +function rand(::Type{HParallelotope}; + N::Type{<:Real}=Float64, + dim::Int=2, + rng::AbstractRNG=GLOBAL_RNG, + seed::Union{Int,Nothing}=nothing) + rng = reseed!(rng, seed) + + while true + D = randn(rng, N, dim, dim) + offset = randn(rng, N, 2 * dim) + + # make sure that the set is not empty: + # `offset[i] >= -offset[i-dim]` for all `i ∈ dim+1:2*dim` + @inbounds for i in (dim + 1):(2 * dim) + offset[i] = -offset[i - dim] + abs(offset[i]) + end + + # convert to polyhedron to check boundedness + clist = _constraints_list_hparallelotope(D, offset, N, typeof(offset)) + Q = HPolyhedron(clist) + if isbounded(Q) + return P = HParallelotope(D, offset; check_consistency=false) + end + # set is unbounded; sample a new set in the next iteration + end +end diff --git a/src/Sets/HParallelotope/volume.jl b/src/Sets/HParallelotope/volume.jl new file mode 100644 index 0000000000..0523ce62bf --- /dev/null +++ b/src/Sets/HParallelotope/volume.jl @@ -0,0 +1,26 @@ +""" + volume(P::HParallelotope) + +Return the volume of a parallelotope in constraint representation. + +### Input + +- `P` -- parallelotope in constraint representation + +### Output + +The volume. + +### Algorithm + +The volume of an ``n``-dimensional parallelotope `P` is ``2^n · |\\det(G)|``, +where ``G`` is the generator matrix of `P`. This can be seen as follows: +The generator matrix transforms the ``n``-dimensional hypercube ``[0, 1]^n`` +to a parallelotope of volume ``|\\det(G)|``. Since the representation of a +parallelotope instead transforms the hypercube ``[-1, 1]^n``, this result has to +be doubled for each dimension. +""" +function volume(P::HParallelotope) + n = dim(P) + return 2^n * abs(det(genmat(P))) +end From 097a2af30124287bf93f6a9b8adf0f81c08dfecf Mon Sep 17 00:00:00 2001 From: schillic Date: Fri, 12 Jul 2024 19:42:02 +0200 Subject: [PATCH 4/6] reorder includes --- src/LazySets.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/LazySets.jl b/src/LazySets.jl index 165abc2712..f07abb1f54 100644 --- a/src/LazySets.jl +++ b/src/LazySets.jl @@ -119,6 +119,11 @@ include("Sets/Ellipsoid/EllipsoidModule.jl") include("Sets/DensePolynomialZonotope/DensePolynomialZonotopeModule.jl") @reexport using ..DensePolynomialZonotopeModule: DensePolynomialZonotope +include("Sets/HPolygon.jl") +include("Sets/HPolygonOpt.jl") +include("Sets/HPolytope.jl") +include("Sets/HPolyhedron.jl") + include("Sets/HParallelotope/HParallelotopeModule.jl") @reexport using ..HParallelotopeModule: HParallelotope, directions, @@ -126,10 +131,6 @@ include("Sets/HParallelotope/HParallelotopeModule.jl") extremal_vertices, offset -include("Sets/HPolygon.jl") -include("Sets/HPolygonOpt.jl") -include("Sets/HPolytope.jl") -include("Sets/HPolyhedron.jl") include("Sets/Hyperplane.jl") include("Sets/Hyperrectangle.jl") From e597be1727d6e8201b1921e57e40c5f8b74e995c Mon Sep 17 00:00:00 2001 From: schillic Date: Fri, 12 Jul 2024 20:43:09 +0200 Subject: [PATCH 5/6] fix docstring --- src/Sets/HParallelotope/HParallelotope.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Sets/HParallelotope/HParallelotope.jl b/src/Sets/HParallelotope/HParallelotope.jl index 55aba9239c..d6e6313c14 100644 --- a/src/Sets/HParallelotope/HParallelotope.jl +++ b/src/Sets/HParallelotope/HParallelotope.jl @@ -20,7 +20,7 @@ form. The `HParallelotope` type represents parallelotopes in constraint form. Let ``D ∈ ℝ^{n × n}`` be a matrix and let ``c ∈ ℝ^{2n}`` be a vector. The parallelotope ``P ⊂ ℝ^n`` generated by the directions matrix ``D`` and the offset vector ``c`` is given by the set of points -``x ∈ ℝ^`` such that: +``x ∈ ℝ^n`` such that: ```math D_i ⋅ x ≤ c_{i},\\text{ and } -D_i ⋅ x ≤ c_{n+i} From 2714f9e55adc1a40098002c45f93d7a58aadf656 Mon Sep 17 00:00:00 2001 From: schillic Date: Fri, 12 Jul 2024 20:51:07 +0200 Subject: [PATCH 6/6] add tests --- src/Sets/HParallelotope/rand.jl | 2 +- test/Sets/HParallelotope.jl | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/Sets/HParallelotope/rand.jl b/src/Sets/HParallelotope/rand.jl index 0a0a6b673f..fc72cfd5d3 100644 --- a/src/Sets/HParallelotope/rand.jl +++ b/src/Sets/HParallelotope/rand.jl @@ -53,5 +53,5 @@ function rand(::Type{HParallelotope}; return P = HParallelotope(D, offset; check_consistency=false) end # set is unbounded; sample a new set in the next iteration - end + end # COV_EXCL_LINE end diff --git a/test/Sets/HParallelotope.jl b/test/Sets/HParallelotope.jl index 95d6035c2f..b579ab5537 100644 --- a/test/Sets/HParallelotope.jl +++ b/test/Sets/HParallelotope.jl @@ -39,6 +39,8 @@ for N in [Float32, Float64, Rational{Int}] HalfSpace(-D[1, :], c[4]), HalfSpace(-D[2, :], c[5]), HalfSpace(-D[3, :], c[6])] + P3 = HParallelotope(zeros(N, 0, 0), zeros(N, 0); check_consistency=false) + @test isempty(constraints_list(P3)) P = HParallelotope(N[1 0; 0 1], N[1, 1, 1, 1]) @@ -79,3 +81,5 @@ for N in [Float64] P = HParallelotope(N[1 0; 0 1], N[1, 1, 1, 1]) @test collect(generators(P)) == [N[1, 0], N[0, 1]] end + +@test !isoperationtype(HParallelotope)