diff --git a/docs/make.jl b/docs/make.jl index 9a368084d2..f382bc7d38 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -116,9 +116,7 @@ makedocs(; sitename="LazySets.jl", "lib/concrete_binary_operations/convex_hull.md", "lib/concrete_binary_operations/difference.md", "lib/concrete_binary_operations/distance.md", - "lib/concrete_binary_operations/exact_sum.md", "lib/concrete_binary_operations/intersection.md", - "lib/concrete_binary_operations/linear_combination.md", "lib/concrete_binary_operations/minkowski_sum.md", "lib/concrete_binary_operations/minkowski_difference.md", "lib/concrete_binary_operations/isdisjoint.md", diff --git a/docs/src/lib/concrete_binary_operations/cartesian_product.md b/docs/src/lib/concrete_binary_operations/cartesian_product.md index f5b7a41359..1ee06de401 100644 --- a/docs/src/lib/concrete_binary_operations/cartesian_product.md +++ b/docs/src/lib/concrete_binary_operations/cartesian_product.md @@ -10,9 +10,6 @@ CurrentModule = LazySets # Cartesian Product ```@docs -cartesian_product(::VPolytope, ::VPolytope) cartesian_product(::LazySet, ::LazySet) -cartesian_product(::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZonotope) -cartesian_product(::SparsePolynomialZonotope, ::SparsePolynomialZonotope) cartesian_product(::SparsePolynomialZonotope, ::AbstractZonotope) ``` diff --git a/docs/src/lib/concrete_binary_operations/convex_hull.md b/docs/src/lib/concrete_binary_operations/convex_hull.md index 074f441aa7..d54242d6ca 100644 --- a/docs/src/lib/concrete_binary_operations/convex_hull.md +++ b/docs/src/lib/concrete_binary_operations/convex_hull.md @@ -12,9 +12,6 @@ CurrentModule = LazySets ```@docs convex_hull(::LazySet, ::LazySet) convex_hull(::HPoly, ::HPoly) -convex_hull(::VPolytope, ::VPolytope) -convex_hull(::VPolygon, ::VPolygon) convex_hull(::Vector{VN}) where {N, VN<:AbstractVector{N}} -convex_hull(::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZonotope) monotone_chain! ``` diff --git a/docs/src/lib/concrete_binary_operations/exact_sum.md b/docs/src/lib/concrete_binary_operations/exact_sum.md deleted file mode 100644 index d193f92c38..0000000000 --- a/docs/src/lib/concrete_binary_operations/exact_sum.md +++ /dev/null @@ -1,14 +0,0 @@ -```@contents -Pages = ["exact_sum.md"] -Depth = 3 -``` - -```@meta -CurrentModule = LazySets -``` - -# Exact Sum - -```@docs -exact_sum(::SparsePolynomialZonotope, ::SparsePolynomialZonotope) -``` diff --git a/docs/src/lib/concrete_binary_operations/intersection.md b/docs/src/lib/concrete_binary_operations/intersection.md index 5f72205673..04afff7a92 100644 --- a/docs/src/lib/concrete_binary_operations/intersection.md +++ b/docs/src/lib/concrete_binary_operations/intersection.md @@ -11,7 +11,6 @@ CurrentModule = LazySets ```@docs intersection(::AbstractSingleton, ::LazySet) -intersection(::Line2D, ::Line2D) intersection(::AbstractHyperrectangle, ::AbstractHyperrectangle) intersection(::Interval, ::HalfSpace) intersection(::Interval, ::Hyperplane) @@ -19,7 +18,6 @@ intersection(::Interval, ::LazySet) intersection(::AbstractHPolygon, ::AbstractHPolygon) intersection(::AbstractPolyhedron{N}, ::AbstractPolyhedron{N}) where {N} intersection(::Union{VPolytope, VPolygon}, ::Union{VPolytope, VPolygon}) -intersection(::VPolygon, ::VPolygon; ::Bool=true) intersection(::UnionSet, ::LazySet) intersection(::UnionSetArray, ::LazySet) intersection(::Universe, ::LazySet) diff --git a/docs/src/lib/concrete_binary_operations/linear_combination.md b/docs/src/lib/concrete_binary_operations/linear_combination.md deleted file mode 100644 index bf7757127f..0000000000 --- a/docs/src/lib/concrete_binary_operations/linear_combination.md +++ /dev/null @@ -1,14 +0,0 @@ -```@contents -Pages = ["linear_combination.md"] -Depth = 3 -``` - -```@meta -CurrentModule = LazySets -``` - -# Linear Combination - -```@docs -linear_combination(::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZonotope) -``` diff --git a/docs/src/lib/concrete_binary_operations/minkowski_sum.md b/docs/src/lib/concrete_binary_operations/minkowski_sum.md index 86af6e4dcd..3ac7970c4f 100644 --- a/docs/src/lib/concrete_binary_operations/minkowski_sum.md +++ b/docs/src/lib/concrete_binary_operations/minkowski_sum.md @@ -12,12 +12,8 @@ CurrentModule = LazySets ```@docs minkowski_sum(::LazySet, ::LazySet) minkowski_sum(::AbstractPolyhedron, ::AbstractPolyhedron) -minkowski_sum(::VPolytope, ::VPolytope) minkowski_sum(::AbstractHyperrectangle, ::AbstractHyperrectangle) minkowski_sum(::AbstractZonotope, ::AbstractZonotope) -minkowski_sum(::VPolygon, ::VPolygon) minkowski_sum(::DensePolynomialZonotope, ::AbstractZonotope) minkowski_sum(::AbstractSingleton, ::AbstractSingleton) -minkowski_sum(::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZonotope) -minkowski_sum(::SparsePolynomialZonotope, ::SparsePolynomialZonotope) ``` diff --git a/docs/src/lib/sets/Line2D.md b/docs/src/lib/sets/Line2D.md index 7c87563a2d..d11bf088f7 100644 --- a/docs/src/lib/sets/Line2D.md +++ b/docs/src/lib/sets/Line2D.md @@ -23,6 +23,7 @@ rand(::Type{Line2D}) project(::AbstractVector, ::Line2D) σ(::AbstractVector, ::Line2D) translate(::Line2D, ::AbstractVector) +intersection(::Line2D, ::Line2D) ``` ```@meta diff --git a/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md b/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md index 40594b0c21..109dfcf6b1 100644 --- a/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md +++ b/docs/src/lib/sets/SimpleSparsePolynomialZonotope.md @@ -21,6 +21,10 @@ remove_redundant_generators(::SimpleSparsePolynomialZonotope) linear_map(::AbstractMatrix, ::SimpleSparsePolynomialZonotope) quadratic_map(::Vector{MT}, ::SimpleSparsePolynomialZonotope) where {N, MT<:AbstractMatrix{N}} quadratic_map(::Vector{MT}, ::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZonotope) where {N, MT<:AbstractMatrix{N}} +cartesian_product(::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZonotope) +convex_hull(::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZonotope) +linear_combination(::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZonotope) +minkowski_sum(::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZonotope) ``` ```@meta diff --git a/docs/src/lib/sets/SparsePolynomialZonotope.md b/docs/src/lib/sets/SparsePolynomialZonotope.md index ec70630afa..a111d668f7 100644 --- a/docs/src/lib/sets/SparsePolynomialZonotope.md +++ b/docs/src/lib/sets/SparsePolynomialZonotope.md @@ -24,6 +24,9 @@ linear_map(::AbstractMatrix, ::SparsePolynomialZonotope) reduce_order(::SparsePolynomialZonotope, ::Real, ::AbstractReductionMethod=GIR05()) ρ(::AbstractVector, ::SparsePolynomialZonotope) translate(::SparsePolynomialZonotope, ::AbstractVector) +cartesian_product(::SparsePolynomialZonotope, ::SparsePolynomialZonotope) +exact_sum(::SparsePolynomialZonotope, ::SparsePolynomialZonotope) +minkowski_sum(::SparsePolynomialZonotope, ::SparsePolynomialZonotope) ``` ```@meta diff --git a/docs/src/lib/sets/VPolygon.md b/docs/src/lib/sets/VPolygon.md index b2c9291e8c..4526834f16 100644 --- a/docs/src/lib/sets/VPolygon.md +++ b/docs/src/lib/sets/VPolygon.md @@ -32,6 +32,9 @@ permute(::VPolygon, ::AbstractVector{Int}) σ(::AbstractVector, ::VPolygon) translate(::VPolygon, ::AbstractVector) translate!(::VPolygon, ::AbstractVector) +convex_hull(::VPolygon, ::VPolygon) +intersection(::VPolygon, ::VPolygon; ::Bool=true) +minkowski_sum(::VPolygon, ::VPolygon) ``` ```@meta diff --git a/docs/src/lib/sets/VPolytope.md b/docs/src/lib/sets/VPolytope.md index 388de5fa3e..0456160783 100644 --- a/docs/src/lib/sets/VPolytope.md +++ b/docs/src/lib/sets/VPolytope.md @@ -32,6 +32,9 @@ linear_map(::AbstractMatrix, ::VPolytope) σ(::AbstractVector, ::VPolytope) translate(::VPolytope, ::AbstractVector) translate!(::VPolytope, ::AbstractVector) +cartesian_product(::VPolytope, ::VPolytope) +convex_hull(::VPolytope, ::VPolytope) +minkowski_sum(::VPolytope, ::VPolytope) ``` ```@meta diff --git a/src/ConcreteOperations/cartesian_product.jl b/src/ConcreteOperations/cartesian_product.jl index 2d05e05a0b..741a2610d7 100644 --- a/src/ConcreteOperations/cartesian_product.jl +++ b/src/ConcreteOperations/cartesian_product.jl @@ -66,50 +66,6 @@ function cartesian_product(X::LazySet, Y::LazySet; backend=nothing, return Pout end -""" - cartesian_product(P1::VPolytope, P2::VPolytope; [backend]=nothing) - -Compute the Cartesian product of two polytopes in vertex representation. - -### Input - -- `P1` -- polytope in vertex representation -- `P2` -- polytope in vertex representation -- `backend` -- (optional, default: `nothing`) backend for polyhedral computation - -### Output - -The `VPolytope` obtained by the concrete Cartesian product of `P1` and `P2`. - -### Notes - -For further information on the supported backends see -[Polyhedra's documentation](https://juliapolyhedra.github.io/). -""" -function cartesian_product(P1::VPolytope, P2::VPolytope; backend=nothing) - require(@__MODULE__, :Polyhedra; fun_name="cartesian_product") - - return _cartesian_product_vrep(P1, P2; backend1=backend, backend2=backend) -end - -function _cartesian_product_vrep(P1, P2; backend1=nothing, backend2=nothing) - if isnothing(backend1) - backend1 = default_polyhedra_backend(P1) - end - if isnothing(backend2) - backend2 = default_polyhedra_backend(P2) - end - - P1′ = polyhedron(P1; backend=backend1) - P2′ = polyhedron(P2; backend=backend2) - Pout = Polyhedra.vcartesianproduct(P1′, P2′) - return VPolytope(Pout) -end - -function cartesian_product(U1::Universe, U2::Universe) - return Universe(dim(U1) + dim(U2)) -end - function _cartesian_product_hrep(X::S1, Y::S2) where {S1<:LazySet,S2<:LazySet} N = promote_type(eltype(X), eltype(Y)) U1 = Universe{N}(dim(X)) @@ -200,67 +156,6 @@ function cartesian_product(U::Universe, H::HalfSpace) return HalfSpace(prepend_zeros(H.a, dim(U)), H.b) end -""" - cartesian_product(P1::SimpleSparsePolynomialZonotope, - P2::SimpleSparsePolynomialZonotope) - -Compute the Cartesian product of two simple sparse polynomial zonotopes. - -### Input - -- `P1` -- simple sparse polynomial zonotope -- `P2` -- simple sparse polynomial zonotope - -### Output - -The Cartesian product of `P1` and `P2`. - -### Algorithm - -This method implements Proposition 3.1.22 in [1]. - -[1] Kochdumper, Niklas. *Extensions of polynomial zonotopes and their application - to verification of cyber-physical systems.* PhD diss., Technische Universität - München, 2022. -""" -function cartesian_product(P1::SimpleSparsePolynomialZonotope, - P2::SimpleSparsePolynomialZonotope) - c = vcat(center(P1), center(P2)) - G = cat(genmat(P1), genmat(P2); dims=(1, 2)) - E = cat(expmat(P1), expmat(P2); dims=(1, 2)) - return SimpleSparsePolynomialZonotope(c, G, E) -end - -""" - cartesian_product(P1::SparsePolynomialZonotope, P2::SparsePolynomialZonotope) - -Compute the Cartesian product of two sparse polynomial zonotopes. - -### Input - -- `P1` -- sparse polynomial zonotope -- `P2` -- sparse polynomial zonotope - -### Output - -The Cartesian product of `P1` and `P2`. - -### Algorithm - -This method implements Proposition 3.1.22 in [1]. - -[1] Kochdumper, Niklas. *Extensions of polynomial zonotopes and their application - to verification of cyber-physical systems.* PhD diss., Technische Universität - München, 2022. -""" -function cartesian_product(P1::SparsePolynomialZonotope, P2::SparsePolynomialZonotope) - c = vcat(center(P1), center(P2)) - G = cat(genmat_dep(P1), genmat_dep(P2); dims=(1, 2)) - GI = cat(genmat_indep(P1), genmat_indep(P2); dims=(1, 2)) - E = cat(expmat(P1), expmat(P2); dims=(1, 2)) - return SparsePolynomialZonotope(c, G, GI, E) -end - """ cartesian_product(SPZ::SparsePolynomialZonotope, Z::AbstractZonotope) diff --git a/src/ConcreteOperations/convex_hull.jl b/src/ConcreteOperations/convex_hull.jl index c10c434f1d..9227b4f1a8 100644 --- a/src/ConcreteOperations/convex_hull.jl +++ b/src/ConcreteOperations/convex_hull.jl @@ -465,119 +465,6 @@ function monotone_chain!(points::Vector{VN}; sort::Bool=true) where {N,VN<:Abstr return resize!(points, m) end -""" - convex_hull(P::VPolygon, Q::VPolygon; [algorithm]::String="monotone_chain") - -Return the convex hull of two polygons in vertex representation. - -### Input - -- `P` -- polygon in vertex representation -- `Q` -- polygon in vertex representation -- `algorithm` -- (optional, default: "monotone_chain") the algorithm used to - compute the convex hull - -### Output - -A new polygon such that its vertices are the convex hull of the two polygons. - -### Notes - -The vertices of the output polygon are sorted in counter-clockwise fashion. -""" -function convex_hull(P::VPolygon, Q::VPolygon; - algorithm::String="monotone_chain") - vunion = [P.vertices; Q.vertices] - convex_hull!(vunion; algorithm=algorithm) - return VPolygon(vunion; apply_convex_hull=false) -end - -""" - convex_hull(P1::VPolytope, P2::VPolytope; [backend]=nothing) - -Compute the convex hull of two polytopes in vertex representation. - -### Input - -- `P1` -- polytope in vertex representation -- `P2` -- polytope in vertex representation -- `backend` -- (optional, default: `nothing`) the polyhedral computation backend - -### Output - -The `VPolytope` obtained by the concrete convex hull of `P1` and `P2`. - -### Notes - -This function takes the union of the vertices of each polytope and then relies -on a concrete convex-hull algorithm. -For low dimensions, a specialized implementation for polygons is used. -For higher dimensions, `convex_hull` relies on the polyhedral backend that can -be specified using the `backend` keyword argument. - -For performance reasons, it is suggested to use the `CDDLib.Library()` backend. -""" -function convex_hull(P1::VPolytope, P2::VPolytope; backend=nothing) - vunion = [P1.vertices; P2.vertices] - convex_hull!(vunion; backend=backend) - return VPolytope(vunion) -end - -""" - convex_hull(P1::HPoly, P2::HPoly; - [backend]=default_polyhedra_backend(P1)) - -Compute the convex hull of the set union of two polyhedra in constraint -representation. - -### Input - -- `P1` -- polyhedron -- `P2` -- polyhedron -- `backend` -- (optional, default: `default_polyhedra_backend(P1)`) the - backend for polyhedral computations - -### Output - -The `HPolyhedron` (resp. `HPolytope`) obtained by the concrete convex hull of -`P1` and `P2`. - -### Notes - -For performance reasons, it is suggested to use the `CDDLib.Library()` backend -for the `convex_hull`. - -For further information on the supported backends see -[Polyhedra's documentation](https://juliapolyhedra.github.io/). -""" -function convex_hull(P1::HPoly, P2::HPoly; - backend=default_polyhedra_backend(P1)) - require(@__MODULE__, :Polyhedra; fun_name="convex_hull") - Pch = Polyhedra.convexhull(polyhedron(P1; backend=backend), - polyhedron(P2; backend=backend)) - removehredundancy!(Pch) - return convert(basetype(P1), Pch) -end - @commutative function convex_hull(X::LazySet, ::EmptySet) return convex_hull(X) end - -""" - convex_hull(P1::SimpleSparsePolynomialZonotope, - P2::SimpleSparsePolynomialZonotope) - -Compute the convex hull of two simple sparse polynomial zonotopes. - -### Input - -- `P1` : simple sparse polynomial zonotopes -- `P2` : simple sparse polynomial zonotopes - -### Output - -Tightest convex simple sparse polynomial zonotope containing `P1` and `P2`. -""" -function convex_hull(P1::SimpleSparsePolynomialZonotope, P2::SimpleSparsePolynomialZonotope) - return linear_combination(linear_combination(P1, P1), linear_combination(P2, P2)) -end diff --git a/src/ConcreteOperations/intersection.jl b/src/ConcreteOperations/intersection.jl index 3f2d8f4bac..accf3c8505 100644 --- a/src/ConcreteOperations/intersection.jl +++ b/src/ConcreteOperations/intersection.jl @@ -43,47 +43,6 @@ function intersection(S1::AbstractSingleton, S2::AbstractSingleton) return _isapprox(element(S1), element(S2)) ? S1 : EmptySet{N}(dim(S1)) end -""" - intersection(L1::Line2D, L2::Line2D) - -Compute the intersection of two two-dimensional lines. - -### Input - -- `L1` -- line -- `L2` -- line - -### Output - -Three outcomes are possible: - -- If the lines are identical, the result is the first line. -- If the lines are parallel and not identical, the result is the empty set. -- Otherwise the result is the set with the unique intersection point. - -### Algorithm - -We first check whether the lines are parallel. -If not, we use [Cramer's rule](https://en.wikipedia.org/wiki/Cramer%27s_rule) -to compute the intersection point. - -### Examples - -The line ``y = x`` intersected with the line ``y = -x + 1`` respectively with -itself: - -```jldoctest -julia> intersection(Line2D([-1.0, 1], 0.0), Line2D([1.0, 1], 1.0)) -Singleton{Float64, Vector{Float64}}([0.5, 0.5]) - -julia> intersection(Line2D([1.0, 1], 1.0), Line2D([1.0, 1], 1.0)) -Line2D{Float64, Vector{Float64}}([1.0, 1.0], 1.0) -``` -""" -function intersection(L1::Line2D, L2::Line2D) - return _intersection_line2d(L1, L2) -end - # this method can also be called with `HalfSpace` arguments function _intersection_line2d(L1, L2) det = right_turn(L1.a, L2.a) @@ -675,42 +634,6 @@ function intersection(P1::Union{VPolygon,VPolytope}, return VPolytope(Pint) end -""" - intersection(P1::VPolygon, P2::VPolygon; apply_convex_hull::Bool=true) - -Compute the intersection of two polygons in vertex representation. - -### Input - -- `P1` -- polygon in vertex representation -- `P2` -- polygon in vertex representation -- `apply_convex_hull` -- (default, optional: `true`) if `false`, skip the - computation of the convex hull of the resulting polygon - -### Output - -A `VPolygon`, or an `EmptySet` if the intersection is empty. - -### Algorithm - -This function applies the [Sutherland–Hodgman polygon clipping -algorithm](https://en.wikipedia.org/wiki/Sutherland%E2%80%93Hodgman_algorithm). -The implementation is based on the one found in -[rosetta code](http://www.rosettacode.org/wiki/Sutherland-Hodgman_polygon_clipping#Julia). -""" -function intersection(P1::VPolygon, P2::VPolygon; apply_convex_hull::Bool=true) - v1 = vertices_list(P1) - v2 = vertices_list(P2) - v12 = _intersection_vrep_2d(v1, v2) - - if isempty(v12) - N = promote_type(eltype(P1), eltype(P2)) - return EmptySet{N}(2) - else - return VPolygon(v12; apply_convex_hull=apply_convex_hull) - end -end - """ intersection(cup::UnionSet, X::LazySet) @@ -830,7 +753,6 @@ The set `X`. end # disambiguations -intersection(U::Universe, ::Universe) = U @commutative intersection(U::Universe, P::AbstractPolyhedron) = P @commutative intersection(U::Universe, S::AbstractSingleton) = S @commutative intersection(U::Universe, X::Interval) = X diff --git a/src/ConcreteOperations/linear_combination.jl b/src/ConcreteOperations/linear_combination.jl index e22edc91e1..09d2875166 100644 --- a/src/ConcreteOperations/linear_combination.jl +++ b/src/ConcreteOperations/linear_combination.jl @@ -1,45 +1,3 @@ -""" - linear_combination(P1::SimpleSparsePolynomialZonotope, - P2::SimpleSparsePolynomialZonotope) - -Compute the linear combination of two simple sparse polynomial zonotopes. - -### Input - -- `P1` -- simple sparse polynomial zonotope -- `P2` -- simple sparse polynomial zonotope - -### Output - -A `SimpleSparsePolynomialZonotope` representing the linear combination of `P1` -and `P2`. - -### Notes - -This method implements the algorithm described in Proposition 3.1.25 of [1]. - -[1] N. Kochdumper. *Extensions of polynomial zonotopes and their application to -verification of cyber-physical systems*. 2021. -""" -function linear_combination(P1::SimpleSparsePolynomialZonotope, - P2::SimpleSparsePolynomialZonotope) - c1, c2 = center(P1), center(P2) - G1, G2 = genmat(P1), genmat(P2) - E1, E2 = expmat(P1), expmat(P2) - - c = (c1 + c2) / 2 - G = hcat(c1 - c2, G1, G1, G2, -G2) / 2 - - N = promote_type(eltype(E1), eltype(E2)) - E = hcat(vcat(zeros(N, nparams(P1) + nparams(P2)), one(N)), - vcat(E1, zeros(N, nparams(P2) + 1, ngens(P1))), - vcat(E1, zeros(N, nparams(P2), ngens(P1)), ones(N, 1, ngens(P1))), - vcat(zeros(N, nparams(P1), ngens(P2)), E2, zeros(N, 1, ngens(P2))), - vcat(zeros(N, nparams(P1), ngens(P2)), E2, ones(N, 1, ngens(P2)))) - - return SimpleSparsePolynomialZonotope(c, G, E) -end - function linear_combination(X::ConvexSet, Y::ConvexSet) return convex_hull(X, Y) end diff --git a/src/ConcreteOperations/minkowski_sum.jl b/src/ConcreteOperations/minkowski_sum.jl index c877157267..332c83ce68 100644 --- a/src/ConcreteOperations/minkowski_sum.jl +++ b/src/ConcreteOperations/minkowski_sum.jl @@ -125,12 +125,6 @@ function minkowski_sum(P::AbstractPolytope, Q::AbstractPolytope; return _minkowski_sum_hrep_preprocess(P, Q, backend, algorithm, prune) end -function minkowski_sum(P::HPolytope, Q::HPolytope; - backend=nothing, algorithm=nothing, prune=true) - res = _minkowski_sum_hrep_preprocess(P, Q, backend, algorithm, prune) - return convert(HPolytope, res) -end - # common code before calling _minkowski_sum_hrep function _minkowski_sum_hrep_preprocess(P, Q, backend, algorithm, prune) require(@__MODULE__, :Polyhedra; fun_name="minkowski_sum") @@ -251,33 +245,6 @@ function minkowski_sum(X::AbstractSingleton, Y::AbstractSingleton) return Singleton(element(X) + element(Y)) end -""" - minkowski_sum(P::VPolygon, Q::VPolygon) - -The Minkowski Sum of two polygons in vertex representation. - -### Input - -- `P` -- polygon in vertex representation -- `Q` -- polygon in vertex representation - -### Output - -A polygon in vertex representation. - -### Algorithm - -We treat each edge of the polygons as a vector, attaching them in polar order -(attaching the tail of the next vector to the head of the previous vector). The -resulting polygonal chain will be a polygon, which is the Minkowski sum of the -given polygons. This algorithm assumes that the vertices of `P` and `Q` are -sorted in counter-clockwise fashion and has linear complexity ``O(m+n)``, where -``m`` and ``n`` are the number of vertices of `P` and `Q`, respectively. -""" -function minkowski_sum(P::VPolygon, Q::VPolygon) - R = _minkowski_sum_vrep_2d(P.vertices, Q.vertices) - return VPolygon(R) -end function _minkowski_sum_vpolygon(P::LazySet, Q::LazySet) return minkowski_sum(convert(VPolygon, P), convert(VPolygon, Q)) @@ -343,46 +310,6 @@ function _minkowski_sum_vrep_2d_singleton(vlistP::Vector{VT}, return R end -""" - minkowski_sum(P1::VPolytope, P2::VPolytope; - [apply_convex_hull]=true, - [backend]=nothing, - [solver]=nothing) - -Compute the Minkowski sum of two polytopes in vertex representation. - -### Input - -- `P1` -- polytope -- `P2` -- polytope -- `apply_convex_hull` -- (optional, default: `true`) if `true`, post-process the - pairwise sums using a convex-hull algorithm -- `backend` -- (optional, default: `nothing`) the backend for - polyhedral computations used to post-process with a - convex hull; see `default_polyhedra_backend(P1)` -- `solver` -- (optional, default: `nothing`) the backend used to - solve the linear program; see - `default_lp_solver_polyhedra(N)` - -### Output - -A new polytope in vertex representation whose vertices are the convex hull of -the sum of all possible sums of vertices of `P1` and `P2`. -""" -function minkowski_sum(P1::VPolytope, P2::VPolytope; - apply_convex_hull::Bool=true, backend=nothing, - solver=nothing) - @assert dim(P1) == dim(P2) "cannot compute the Minkowski sum of two " * - "polytopes of dimension $(dim(P1)) and $(dim(P2)), respectively" - - vlist1 = _vertices_list(P1, backend) - vlist2 = _vertices_list(P2, backend) - Vout = _minkowski_sum_vrep_nd(vlist1, vlist2; - apply_convex_hull=apply_convex_hull, - backend=backend, solver=solver) - return VPolytope(Vout) -end - function _minkowski_sum_vrep_nd(vlist1::Vector{VT}, vlist2::Vector{VT}; apply_convex_hull::Bool=true, backend=nothing, solver=nothing) where {N,VT<:AbstractVector{N}} @@ -436,59 +363,6 @@ for T in [:LazySet, :AbstractPolyhedron, :AbstractPolytope, :AbstractZonotope, end end -""" - minkowski_sum(P1::SimpleSparsePolynomialZonotope, - P2::SimpleSparsePolynomialZonotope) - -Compute the Minkowski sum of two simple sparse polynomial zonotopes. - -### Input - -- `P1` -- simple sparse polynomial zonotope -- `P2` -- simple sparse polynomial zonotope - -### Output - -The Minkowski sum of `P1` and `P2`. -""" -function minkowski_sum(P1::SimpleSparsePolynomialZonotope, - P2::SimpleSparsePolynomialZonotope) - c = center(P1) + center(P2) - G = hcat(genmat(P1), genmat(P2)) - E = cat(expmat(P1), expmat(P2); dims=(1, 2)) - return SimpleSparsePolynomialZonotope(c, G, E) -end - -""" - minkowski_sum(P1::SparsePolynomialZonotope, P2::SparsePolynomialZonotope) - -Compute the Minkowski sum of two sparse polyomial zonotopes. - -### Input - -- `P1` -- sparse polynomial zonotope -- `P2` -- sparse polynomial zonotope - -### Output - -The Minkowski sum of `P1` and `P2`. - -### Algorithm - -See Proposition 3.1.19 in [1]. - -[1] Kochdumper. *Extensions of polynomial zonotopes and their application to - verification of cyber-physical systems.* PhD diss., TU Munich, 2022. -""" -function minkowski_sum(P1::SparsePolynomialZonotope, - P2::SparsePolynomialZonotope) - c = center(P1) + center(P2) - G = hcat(genmat_dep(P1), genmat_dep(P2)) - GI = hcat(genmat_indep(P1), genmat_indep(P2)) - E = cat(expmat(P1), expmat(P2); dims=(1, 2)) - return SparsePolynomialZonotope(c, G, GI, E) -end - # See Proposition 3.1.19 in [1]. # [1] Kochdumper. *Extensions of polynomial zonotopes and their application to # verification of cyber-physical systems.* PhD diss., TU Munich, 2022. diff --git a/src/LazySets.jl b/src/LazySets.jl index 59e3588aab..dd3233a5fb 100644 --- a/src/LazySets.jl +++ b/src/LazySets.jl @@ -268,7 +268,6 @@ include("ConcreteOperations/cartesian_product.jl") include("ConcreteOperations/convex_hull.jl") include("ConcreteOperations/difference.jl") include("ConcreteOperations/distance.jl") -include("ConcreteOperations/exact_sum.jl") include("ConcreteOperations/intersection.jl") include("ConcreteOperations/isapprox.jl") include("ConcreteOperations/isdisjoint.jl") diff --git a/src/Sets/Ball1/Ball1Module.jl b/src/Sets/Ball1/Ball1Module.jl index 33f7cc8ce0..d55dece5d0 100644 --- a/src/Sets/Ball1/Ball1Module.jl +++ b/src/Sets/Ball1/Ball1Module.jl @@ -20,7 +20,6 @@ export Ball1 include("Ball1.jl") -include("ball_norm.jl") include("center.jl") include("constraints_list.jl") include("high.jl") @@ -28,7 +27,6 @@ include("in.jl") include("isoperationtype.jl") include("low.jl") include("project.jl") -include("radius_ball.jl") include("rand.jl") include("reflect.jl") include("scale.jl") @@ -37,6 +35,9 @@ include("support_vector.jl") include("translate.jl") include("vertices_list.jl") +include("ball_norm.jl") +include("radius_ball.jl") + include("init.jl") end # module diff --git a/src/Sets/Ball2/Ball2Module.jl b/src/Sets/Ball2/Ball2Module.jl index 6ec966e8f0..52ed99d96b 100644 --- a/src/Sets/Ball2/Ball2Module.jl +++ b/src/Sets/Ball2/Ball2Module.jl @@ -20,13 +20,10 @@ export Ball2 include("Ball2.jl") include("area.jl") -include("ball_norm.jl") include("center.jl") -include("chebyshev_center_radius.jl") include("in.jl") include("isoperationtype.jl") include("project.jl") -include("radius_ball.jl") include("rand.jl") include("reflect.jl") include("sample.jl") @@ -38,6 +35,10 @@ include("volume.jl") include("isdisjoint.jl") include("issubset.jl") +include("ball_norm.jl") +include("chebyshev_center_radius.jl") +include("radius_ball.jl") + function ○(c::VN, r::N) where {N<:AbstractFloat,VN<:AbstractVector{N}} return Ball2(c, r) end diff --git a/src/Sets/BallInf/BallInfModule.jl b/src/Sets/BallInf/BallInfModule.jl index 927f61e4a4..19dc67cfe6 100644 --- a/src/Sets/BallInf/BallInfModule.jl +++ b/src/Sets/BallInf/BallInfModule.jl @@ -21,17 +21,11 @@ export BallInf include("BallInf.jl") include("area.jl") -include("ball_norm.jl") include("center.jl") -include("genmat.jl") include("high.jl") -include("isflat.jl") include("isoperationtype.jl") include("low.jl") -include("ngens.jl") include("radius.jl") -include("radius_ball.jl") -include("radius_hyperrectangle.jl") include("rand.jl") include("reflect.jl") include("volume.jl") @@ -41,6 +35,13 @@ include("support_function.jl") include("support_vector.jl") include("translate.jl") +include("ball_norm.jl") +include("genmat.jl") +include("isflat.jl") +include("ngens.jl") +include("radius_ball.jl") +include("radius_hyperrectangle.jl") + include("init.jl") function □(c::VN, r::N) where {N,VN<:AbstractVector{N}} diff --git a/src/Sets/Ballp/BallpModule.jl b/src/Sets/Ballp/BallpModule.jl index a57dbd922d..339ae5d1fe 100644 --- a/src/Sets/Ballp/BallpModule.jl +++ b/src/Sets/Ballp/BallpModule.jl @@ -16,16 +16,17 @@ export Ballp include("Ballp.jl") -include("ball_norm.jl") include("center.jl") include("isoperationtype.jl") include("project.jl") -include("radius_ball.jl") include("rand.jl") include("reflect.jl") include("scale.jl") include("translate.jl") +include("ball_norm.jl") +include("radius_ball.jl") + include("init.jl") end # module diff --git a/src/Sets/Ellipsoid/EllipsoidModule.jl b/src/Sets/Ellipsoid/EllipsoidModule.jl index 5b1b01ef98..49d03ea4a1 100644 --- a/src/Sets/Ellipsoid/EllipsoidModule.jl +++ b/src/Sets/Ellipsoid/EllipsoidModule.jl @@ -24,11 +24,12 @@ include("isoperationtype.jl") include("rand.jl") include("affine_map.jl") include("linear_map.jl") -include("shape_matrix.jl") include("support_function.jl") include("support_vector.jl") include("translate.jl") +include("shape_matrix.jl") + function ○(c::VN, shape_matrix::MN) where {N<:AbstractFloat, VN<:AbstractVector{N}, diff --git a/src/Sets/EmptySet/EmptySetModule.jl b/src/Sets/EmptySet/EmptySetModule.jl index cddde7837b..0e7bf906be 100644 --- a/src/Sets/EmptySet/EmptySetModule.jl +++ b/src/Sets/EmptySet/EmptySetModule.jl @@ -23,39 +23,39 @@ include("EmptySet.jl") include("an_element.jl") include("area.jl") -include("chebyshev_center_radius.jl") include("complement.jl") include("diameter.jl") include("dim.jl") include("high.jl") -include("in.jl") include("isbounded.jl") include("isboundedtype.jl") include("isconvextype.jl") include("isempty.jl") include("isoperationtype.jl") include("isuniversal.jl") -include("linear_map.jl") include("low.jl") include("norm.jl") -include("project.jl") include("radius.jl") include("rand.jl") include("rectify.jl") include("reflect.jl") +include("vertices_list.jl") +include("vertices.jl") +include("volume.jl") +include("in.jl") +include("linear_map.jl") +include("project.jl") include("scale.jl") include("support_function.jl") include("support_vector.jl") include("translate.jl") -include("vertices_list.jl") -include("vertices.jl") -include("volume.jl") - include("convex_hull.jl") include("intersection.jl") include("isdisjoint.jl") include("issubset.jl") +include("chebyshev_center_radius.jl") + """ plot_recipe(∅::EmptySet{N}, [ε]=zero(N)) where {N} diff --git a/src/Sets/HParallelotope/HParallelotopeModule.jl b/src/Sets/HParallelotope/HParallelotopeModule.jl index 19d47f8c39..983ca6983c 100644 --- a/src/Sets/HParallelotope/HParallelotopeModule.jl +++ b/src/Sets/HParallelotope/HParallelotopeModule.jl @@ -24,18 +24,19 @@ export HParallelotope, include("HParallelotope.jl") -include("base_vertex.jl") include("center.jl") include("constraints_list.jl") include("dim.jl") +include("isoperationtype.jl") +include("rand.jl") +include("volume.jl") + +include("base_vertex.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") include("convert.jl") diff --git a/src/Sets/HPolyhedron/HPolyhedronModule.jl b/src/Sets/HPolyhedron/HPolyhedronModule.jl index cbb2ced2ba..49b7846274 100644 --- a/src/Sets/HPolyhedron/HPolyhedronModule.jl +++ b/src/Sets/HPolyhedron/HPolyhedronModule.jl @@ -17,7 +17,7 @@ using ReachabilityBase.Distribution: reseed! using ReachabilityBase.Require: require @reexport import ..API: constraints_list, dim, isempty, isoperationtype, rand, - permute, ρ, σ, translate + permute, ρ, σ, translate, convex_hull @reexport import ..LazySets: is_hyperplanar, normalize, remove_redundant_constraints, remove_redundant_constraints!, tohrep, tovrep, @@ -41,6 +41,7 @@ include("permute.jl") include("support_function.jl") include("support_vector.jl") include("translate.jl") +include("convex_hull.jl") include("is_hyperplanar.jl") include("normalize.jl") diff --git a/src/Sets/HPolyhedron/convex_hull.jl b/src/Sets/HPolyhedron/convex_hull.jl new file mode 100644 index 0000000000..88dfce7f89 --- /dev/null +++ b/src/Sets/HPolyhedron/convex_hull.jl @@ -0,0 +1,35 @@ +""" + convex_hull(P1::HPoly, P2::HPoly; + [backend]=default_polyhedra_backend(P1)) + +Compute the convex hull of the set union of two polyhedra in constraint +representation. + +### Input + +- `P1` -- polyhedron +- `P2` -- polyhedron +- `backend` -- (optional, default: `default_polyhedra_backend(P1)`) the + backend for polyhedral computations + +### Output + +The `HPolyhedron` (resp. `HPolytope`) obtained by the concrete convex hull of +`P1` and `P2`. + +### Notes + +For performance reasons, it is suggested to use the `CDDLib.Library()` backend +for the `convex_hull`. + +For further information on the supported backends see +[Polyhedra's documentation](https://juliapolyhedra.github.io/). +""" +function convex_hull(P1::HPoly, P2::HPoly; + backend=default_polyhedra_backend(P1)) + require(@__MODULE__, :Polyhedra; fun_name="convex_hull") + Pch = Polyhedra.convexhull(LazySets.polyhedron(P1; backend=backend), + LazySets.polyhedron(P2; backend=backend)) + Polyhedra.removehredundancy!(Pch) + return convert(basetype(P1), Pch) +end diff --git a/src/Sets/HPolytope/HPolytopeModule.jl b/src/Sets/HPolytope/HPolytopeModule.jl index f8d16be8a7..5b7c11f3fd 100644 --- a/src/Sets/HPolytope/HPolytopeModule.jl +++ b/src/Sets/HPolytope/HPolytopeModule.jl @@ -4,14 +4,15 @@ using Reexport, Requires using ..LazySets: AbstractPolytope, LazySet, AbstractLinearMapAlgorithm, default_polyhedra_backend, vertices_list_1d, _linear_map_hrep, - _normal_Vector + _minkowski_sum_hrep_preprocess, _normal_Vector using ..HalfSpaceModule: HalfSpace using Random: AbstractRNG, GLOBAL_RNG using ReachabilityBase.Distribution: reseed! using ReachabilityBase.Comparison: isapproxzero, _ztol using ReachabilityBase.Require: require -@reexport import ..API: isbounded, isoperationtype, rand, vertices_list +@reexport import ..API: isbounded, isoperationtype, rand, vertices_list, + minkowski_sum @reexport import ..LazySets: _linear_map_hrep_helper, _vertices_list import Base: convert @reexport using ..API @@ -25,6 +26,7 @@ include("isoperationtype.jl") include("rand.jl") include("vertices_list.jl") include("linear_map.jl") +include("minkowski_sum.jl") include("convert.jl") diff --git a/src/Sets/HPolytope/minkowski_sum.jl b/src/Sets/HPolytope/minkowski_sum.jl new file mode 100644 index 0000000000..a503cc20e7 --- /dev/null +++ b/src/Sets/HPolytope/minkowski_sum.jl @@ -0,0 +1,5 @@ +function minkowski_sum(P::HPolytope, Q::HPolytope; + backend=nothing, algorithm=nothing, prune=true) + res = _minkowski_sum_hrep_preprocess(P, Q, backend, algorithm, prune) + return convert(HPolytope, res) +end diff --git a/src/Sets/Interval/IntervalModule.jl b/src/Sets/Interval/IntervalModule.jl index 24a938569d..3facf15957 100644 --- a/src/Sets/Interval/IntervalModule.jl +++ b/src/Sets/Interval/IntervalModule.jl @@ -27,41 +27,35 @@ export Interval include("Interval.jl") -include("affine_map.jl") include("an_element.jl") include("center.jl") -include("chebyshev_center_radius.jl") include("complement.jl") include("constraints_list.jl") -include("convert.jl") include("convex_hull.jl") include("diameter.jl") include("dim.jl") -include("exponential_map.jl") include("extrema.jl") include("high.jl") -include("in.jl") include("isflat.jl") include("isoperationtype.jl") -include("linear_map.jl") include("low.jl") -include("ngens.jl") include("norm.jl") -include("permute.jl") -include("project.jl") include("radius.jl") -include("radius_hyperrectangle.jl") include("rand.jl") include("rectify.jl") include("reflect.jl") +include("vertices_list.jl") +include("volume.jl") +include("affine_map.jl") +include("exponential_map.jl") +include("in.jl") +include("linear_map.jl") +include("permute.jl") +include("project.jl") include("scale.jl") -include("split.jl") include("support_vector.jl") include("support_function.jl") include("translate.jl") -include("vertices_list.jl") -include("volume.jl") - include("difference.jl") include("distance.jl") include("intersection.jl") @@ -73,6 +67,13 @@ include("issubset.jl") include("minkowski_difference.jl") include("minkowski_sum.jl") +include("chebyshev_center_radius.jl") +include("ngens.jl") +include("radius_hyperrectangle.jl") +include("split.jl") + +include("convert.jl") + """ -(x::Interval, y::Interval) diff --git a/src/Sets/Line2D/Line2DModule.jl b/src/Sets/Line2D/Line2DModule.jl index 78a1c33ded..d4a9610acc 100644 --- a/src/Sets/Line2D/Line2DModule.jl +++ b/src/Sets/Line2D/Line2DModule.jl @@ -3,8 +3,9 @@ module Line2DModule using Reexport, Requires using ..LazySets: AbstractPolyhedron, AbstractLinearMapAlgorithm, - _constraints_list_hyperplane, _linear_map_hrep, - _non_element_halfspace, _σ_hyperplane_halfspace + _constraints_list_hyperplane, _intersection_line2d, + _linear_map_hrep, _non_element_halfspace, + _σ_hyperplane_halfspace using LinearAlgebra: dot using Random: AbstractRNG, GLOBAL_RNG using ReachabilityBase.Arrays: nonzero_indices @@ -14,7 +15,7 @@ using ReachabilityBase.Require: require @reexport import ..API: an_element, constraints_list, dim, isbounded, isempty, isoperationtype, isuniversal, rand, ∈, project, σ, - translate + translate, intersection @reexport import ..LazySets: constrained_dimensions, _linear_map_hrep_helper @reexport using ..API @@ -23,7 +24,6 @@ export Line2D include("Line2D.jl") include("an_element.jl") -include("constrained_dimensions.jl") include("constraints_list.jl") include("dim.jl") include("isbounded.jl") @@ -36,6 +36,9 @@ include("linear_map.jl") include("project.jl") include("support_vector.jl") include("translate.jl") +include("intersection.jl") + +include("constrained_dimensions.jl") include("init.jl") diff --git a/src/Sets/Line2D/intersection.jl b/src/Sets/Line2D/intersection.jl new file mode 100644 index 0000000000..3661416fde --- /dev/null +++ b/src/Sets/Line2D/intersection.jl @@ -0,0 +1,40 @@ +""" + intersection(L1::Line2D, L2::Line2D) + +Compute the intersection of two two-dimensional lines. + +### Input + +- `L1` -- line +- `L2` -- line + +### Output + +Three outcomes are possible: + +- If the lines are identical, the result is the first line. +- If the lines are parallel and not identical, the result is the empty set. +- Otherwise the result is the set with the unique intersection point. + +### Algorithm + +We first check whether the lines are parallel. +If not, we use [Cramer's rule](https://en.wikipedia.org/wiki/Cramer%27s_rule) +to compute the intersection point. + +### Examples + +The line ``y = x`` intersected with the line ``y = -x + 1`` respectively with +itself: + +```jldoctest +julia> intersection(Line2D([-1.0, 1], 0.0), Line2D([1.0, 1], 1.0)) +Singleton{Float64, Vector{Float64}}([0.5, 0.5]) + +julia> intersection(Line2D([1.0, 1], 1.0), Line2D([1.0, 1], 1.0)) +Line2D{Float64, Vector{Float64}}([1.0, 1.0], 1.0) +``` +""" +function intersection(L1::Line2D, L2::Line2D) + return _intersection_line2d(L1, L2) +end diff --git a/src/Sets/LineSegment/LineSegmentModule.jl b/src/Sets/LineSegment/LineSegmentModule.jl index 9203c5fc89..3d0730d153 100644 --- a/src/Sets/LineSegment/LineSegmentModule.jl +++ b/src/Sets/LineSegment/LineSegmentModule.jl @@ -25,10 +25,7 @@ include("an_element.jl") include("center.jl") include("constraints_list.jl") include("dim.jl") -include("generators.jl") -include("genmat.jl") include("isoperationtype.jl") -include("ngens.jl") include("rand.jl") include("vertices_list.jl") include("in.jl") @@ -36,11 +33,13 @@ include("scale.jl") include("support_function.jl") include("support_vector.jl") include("translate.jl") - include("intersection.jl") +include("generators.jl") +include("genmat.jl") include("halfspace_left.jl") include("halfspace_right.jl") +include("ngens.jl") include("init.jl") diff --git a/src/Sets/SimpleSparsePolynomialZonotope/SimpleSparsePolynomialZonotopeModule.jl b/src/Sets/SimpleSparsePolynomialZonotope/SimpleSparsePolynomialZonotopeModule.jl index 947b709af1..a6248ef82a 100644 --- a/src/Sets/SimpleSparsePolynomialZonotope/SimpleSparsePolynomialZonotopeModule.jl +++ b/src/Sets/SimpleSparsePolynomialZonotope/SimpleSparsePolynomialZonotopeModule.jl @@ -9,7 +9,8 @@ using Random: AbstractRNG, GLOBAL_RNG using ReachabilityBase.Distribution: reseed! using ReachabilityBase.Comparison: isapproxzero -@reexport import ..API: convex_hull, center, isoperationtype, rand, linear_map +@reexport import ..API: convex_hull, center, isoperationtype, rand, linear_map, + cartesian_product, linear_combination, minkowski_sum @reexport import ..LazySets: expmat, genmat, genmat_dep, genmat_indep, ngens, ngens_indep, polynomial_order, remove_redundant_generators @@ -26,17 +27,21 @@ const SSPZ = SimpleSparsePolynomialZonotope include("center.jl") include("convex_hull.jl") +include("isoperationtype.jl") +include("rand.jl") +include("linear_map.jl") +include("cartesian_product.jl") +include("linear_combination.jl") +include("minkowski_sum.jl") + include("genmat.jl") include("genmat_dep.jl") include("genmat_indep.jl") -include("isoperationtype.jl") include("ngens.jl") include("ngens_indep.jl") include("polynomial_order.jl") include("remove_redundant_generators.jl") include("expmat.jl") include("quadratic_map.jl") -include("rand.jl") -include("linear_map.jl") end # module diff --git a/src/Sets/SimpleSparsePolynomialZonotope/cartesian_product.jl b/src/Sets/SimpleSparsePolynomialZonotope/cartesian_product.jl new file mode 100644 index 0000000000..b22a60c79b --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/cartesian_product.jl @@ -0,0 +1,30 @@ +""" + cartesian_product(P1::SimpleSparsePolynomialZonotope, + P2::SimpleSparsePolynomialZonotope) + +Compute the Cartesian product of two simple sparse polynomial zonotopes. + +### Input + +- `P1` -- simple sparse polynomial zonotope +- `P2` -- simple sparse polynomial zonotope + +### Output + +The Cartesian product of `P1` and `P2`. + +### Algorithm + +This method implements Proposition 3.1.22 in [1]. + +[1] Kochdumper, Niklas. *Extensions of polynomial zonotopes and their application + to verification of cyber-physical systems.* PhD diss., Technische Universität + München, 2022. +""" +function cartesian_product(P1::SimpleSparsePolynomialZonotope, + P2::SimpleSparsePolynomialZonotope) + c = vcat(center(P1), center(P2)) + G = cat(genmat(P1), genmat(P2); dims=(1, 2)) + E = cat(expmat(P1), expmat(P2); dims=(1, 2)) + return SimpleSparsePolynomialZonotope(c, G, E) +end diff --git a/src/Sets/SimpleSparsePolynomialZonotope/convex_hull.jl b/src/Sets/SimpleSparsePolynomialZonotope/convex_hull.jl index 3914715e42..e41b809db1 100644 --- a/src/Sets/SimpleSparsePolynomialZonotope/convex_hull.jl +++ b/src/Sets/SimpleSparsePolynomialZonotope/convex_hull.jl @@ -14,3 +14,22 @@ The tightest convex simple sparse polynomial zonotope containing `P`. function convex_hull(P::SimpleSparsePolynomialZonotope) return linear_combination(P, P) end + +""" + convex_hull(P1::SimpleSparsePolynomialZonotope, + P2::SimpleSparsePolynomialZonotope) + +Compute the convex hull of two simple sparse polynomial zonotopes. + +### Input + +- `P1` : simple sparse polynomial zonotopes +- `P2` : simple sparse polynomial zonotopes + +### Output + +Tightest convex simple sparse polynomial zonotope containing `P1` and `P2`. +""" +function convex_hull(P1::SimpleSparsePolynomialZonotope, P2::SimpleSparsePolynomialZonotope) + return linear_combination(linear_combination(P1, P1), linear_combination(P2, P2)) +end diff --git a/src/Sets/SimpleSparsePolynomialZonotope/linear_combination.jl b/src/Sets/SimpleSparsePolynomialZonotope/linear_combination.jl new file mode 100644 index 0000000000..762cb581cf --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/linear_combination.jl @@ -0,0 +1,41 @@ +""" + linear_combination(P1::SimpleSparsePolynomialZonotope, + P2::SimpleSparsePolynomialZonotope) + +Compute the linear combination of two simple sparse polynomial zonotopes. + +### Input + +- `P1` -- simple sparse polynomial zonotope +- `P2` -- simple sparse polynomial zonotope + +### Output + +A `SimpleSparsePolynomialZonotope` representing the linear combination of `P1` +and `P2`. + +### Notes + +This method implements the algorithm described in Proposition 3.1.25 of [1]. + +[1] N. Kochdumper. *Extensions of polynomial zonotopes and their application to +verification of cyber-physical systems*. 2021. +""" +function linear_combination(P1::SimpleSparsePolynomialZonotope, + P2::SimpleSparsePolynomialZonotope) + c1, c2 = center(P1), center(P2) + G1, G2 = genmat(P1), genmat(P2) + E1, E2 = expmat(P1), expmat(P2) + + c = (c1 + c2) / 2 + G = hcat(c1 - c2, G1, G1, G2, -G2) / 2 + + N = promote_type(eltype(E1), eltype(E2)) + E = hcat(vcat(zeros(N, nparams(P1) + nparams(P2)), one(N)), + vcat(E1, zeros(N, nparams(P2) + 1, ngens(P1))), + vcat(E1, zeros(N, nparams(P2), ngens(P1)), ones(N, 1, ngens(P1))), + vcat(zeros(N, nparams(P1), ngens(P2)), E2, zeros(N, 1, ngens(P2))), + vcat(zeros(N, nparams(P1), ngens(P2)), E2, ones(N, 1, ngens(P2)))) + + return SimpleSparsePolynomialZonotope(c, G, E) +end diff --git a/src/Sets/SimpleSparsePolynomialZonotope/minkowski_sum.jl b/src/Sets/SimpleSparsePolynomialZonotope/minkowski_sum.jl new file mode 100644 index 0000000000..e9c9c489e2 --- /dev/null +++ b/src/Sets/SimpleSparsePolynomialZonotope/minkowski_sum.jl @@ -0,0 +1,22 @@ +""" + minkowski_sum(P1::SimpleSparsePolynomialZonotope, + P2::SimpleSparsePolynomialZonotope) + +Compute the Minkowski sum of two simple sparse polynomial zonotopes. + +### Input + +- `P1` -- simple sparse polynomial zonotope +- `P2` -- simple sparse polynomial zonotope + +### Output + +The Minkowski sum of `P1` and `P2`. +""" +function minkowski_sum(P1::SimpleSparsePolynomialZonotope, + P2::SimpleSparsePolynomialZonotope) + c = center(P1) + center(P2) + G = hcat(genmat(P1), genmat(P2)) + E = cat(expmat(P1), expmat(P2); dims=(1, 2)) + return SimpleSparsePolynomialZonotope(c, G, E) +end diff --git a/src/Sets/Singleton/SingletonModule.jl b/src/Sets/Singleton/SingletonModule.jl index 5ab75471ec..e2ca545c4b 100644 --- a/src/Sets/Singleton/SingletonModule.jl +++ b/src/Sets/Singleton/SingletonModule.jl @@ -18,7 +18,6 @@ include("Singleton.jl") include("isoperationtype.jl") include("rand.jl") include("rectify.jl") -include("singleton_list.jl") include("linear_map.jl") include("permute.jl") include("project.jl") @@ -26,5 +25,6 @@ include("scale.jl") include("translate.jl") include("element.jl") +include("singleton_list.jl") end # module diff --git a/src/Sets/SparsePolynomialZonotope/SparsePolynomialZonotopeModule.jl b/src/Sets/SparsePolynomialZonotope/SparsePolynomialZonotopeModule.jl index e3599ac82a..1738b6764c 100644 --- a/src/Sets/SparsePolynomialZonotope/SparsePolynomialZonotopeModule.jl +++ b/src/Sets/SparsePolynomialZonotope/SparsePolynomialZonotopeModule.jl @@ -12,7 +12,7 @@ using ReachabilityBase.Distribution: reseed! using ReachabilityBase.Require: require @reexport import ..API: center, extrema, isoperationtype, rand, linear_map, ρ, - translate + translate, cartesian_product, exact_sum, minkowski_sum @reexport import ..LazySets: expmat, genmat_dep, genmat_indep, ngens_dep, ngens_indep, nparams, polynomial_order, reduce_order, remove_redundant_generators @@ -27,20 +27,23 @@ include("SparsePolynomialZonotope.jl") const SPZ = SparsePolynomialZonotope include("center.jl") -include("expmat.jl") include("extrema.jl") -include("genmat_dep.jl") -include("genmat_indep.jl") include("isoperationtype.jl") -include("polynomial_order.jl") include("rand.jl") -include("remove_redundant_generators.jl") include("linear_map.jl") -include("reduce_order.jl") include("support_function.jl") include("translate.jl") +include("cartesian_product.jl") +include("exact_sum.jl") +include("minkowski_sum.jl") +include("expmat.jl") +include("genmat_dep.jl") +include("genmat_indep.jl") include("indexvector.jl") +include("polynomial_order.jl") +include("remove_redundant_generators.jl") +include("reduce_order.jl") include("init.jl") diff --git a/src/Sets/SparsePolynomialZonotope/cartesian_product.jl b/src/Sets/SparsePolynomialZonotope/cartesian_product.jl new file mode 100644 index 0000000000..070e3cb6ac --- /dev/null +++ b/src/Sets/SparsePolynomialZonotope/cartesian_product.jl @@ -0,0 +1,29 @@ +""" + cartesian_product(P1::SparsePolynomialZonotope, P2::SparsePolynomialZonotope) + +Compute the Cartesian product of two sparse polynomial zonotopes. + +### Input + +- `P1` -- sparse polynomial zonotope +- `P2` -- sparse polynomial zonotope + +### Output + +The Cartesian product of `P1` and `P2`. + +### Algorithm + +This method implements Proposition 3.1.22 in [1]. + +[1] Kochdumper, Niklas. *Extensions of polynomial zonotopes and their application + to verification of cyber-physical systems.* PhD diss., Technische Universität + München, 2022. +""" +function cartesian_product(P1::SparsePolynomialZonotope, P2::SparsePolynomialZonotope) + c = vcat(center(P1), center(P2)) + G = cat(genmat_dep(P1), genmat_dep(P2); dims=(1, 2)) + GI = cat(genmat_indep(P1), genmat_indep(P2); dims=(1, 2)) + E = cat(expmat(P1), expmat(P2); dims=(1, 2)) + return SparsePolynomialZonotope(c, G, GI, E) +end diff --git a/src/ConcreteOperations/exact_sum.jl b/src/Sets/SparsePolynomialZonotope/exact_sum.jl similarity index 100% rename from src/ConcreteOperations/exact_sum.jl rename to src/Sets/SparsePolynomialZonotope/exact_sum.jl diff --git a/src/Sets/SparsePolynomialZonotope/minkowski_sum.jl b/src/Sets/SparsePolynomialZonotope/minkowski_sum.jl new file mode 100644 index 0000000000..07606f52bd --- /dev/null +++ b/src/Sets/SparsePolynomialZonotope/minkowski_sum.jl @@ -0,0 +1,29 @@ +""" + minkowski_sum(P1::SparsePolynomialZonotope, P2::SparsePolynomialZonotope) + +Compute the Minkowski sum of two sparse polyomial zonotopes. + +### Input + +- `P1` -- sparse polynomial zonotope +- `P2` -- sparse polynomial zonotope + +### Output + +The Minkowski sum of `P1` and `P2`. + +### Algorithm + +See Proposition 3.1.19 in [1]. + +[1] Kochdumper. *Extensions of polynomial zonotopes and their application to + verification of cyber-physical systems.* PhD diss., TU Munich, 2022. +""" +function minkowski_sum(P1::SparsePolynomialZonotope, + P2::SparsePolynomialZonotope) + c = center(P1) + center(P2) + G = hcat(genmat_dep(P1), genmat_dep(P2)) + GI = hcat(genmat_indep(P1), genmat_indep(P2)) + E = cat(expmat(P1), expmat(P2); dims=(1, 2)) + return SparsePolynomialZonotope(c, G, GI, E) +end diff --git a/src/Sets/Universe/UniverseModule.jl b/src/Sets/Universe/UniverseModule.jl index 7acdeb121e..04a5a0f26e 100644 --- a/src/Sets/Universe/UniverseModule.jl +++ b/src/Sets/Universe/UniverseModule.jl @@ -12,7 +12,7 @@ using ReachabilityBase.Require: require diameter, dim, isbounded, isboundedtype, isempty, isoperationtype, isuniversal, norm, radius, rand, reflect, ∈, permute, project, scale, scale!, ρ, σ, - translate, translate! + translate, translate!, cartesian_product, intersection @reexport import ..LazySets: constrained_dimensions, linear_map_inverse, tosimplehrep @reexport using ..API @@ -43,6 +43,8 @@ include("scale.jl") include("support_function.jl") include("support_vector.jl") include("translate.jl") +include("cartesian_product.jl") +include("intersection.jl") include("constrained_dimensions.jl") include("linear_map_inverse.jl") diff --git a/src/Sets/Universe/cartesian_product.jl b/src/Sets/Universe/cartesian_product.jl new file mode 100644 index 0000000000..8aea30ef0c --- /dev/null +++ b/src/Sets/Universe/cartesian_product.jl @@ -0,0 +1,3 @@ +function cartesian_product(U1::Universe, U2::Universe) + return Universe(dim(U1) + dim(U2)) +end diff --git a/src/Sets/Universe/intersection.jl b/src/Sets/Universe/intersection.jl new file mode 100644 index 0000000000..778b34f393 --- /dev/null +++ b/src/Sets/Universe/intersection.jl @@ -0,0 +1 @@ +intersection(U::Universe, ::Universe) = U diff --git a/src/Sets/VPolygon/VPolygonModule.jl b/src/Sets/VPolygon/VPolygonModule.jl index 7750eb0d13..4ac8645ebc 100644 --- a/src/Sets/VPolygon/VPolygonModule.jl +++ b/src/Sets/VPolygon/VPolygonModule.jl @@ -2,8 +2,9 @@ module VPolygonModule using Reexport, Requires -using ..LazySets: AbstractPolygon, LazySet, AbstractHPolygon, convex_hull, - halfspace_left, is_right_turn, _area_vlist, _linear_map_vrep +using ..LazySets: AbstractPolygon, LazySet, AbstractHPolygon, halfspace_left, + is_right_turn, _area_vlist, _intersection_vrep_2d, + _linear_map_vrep, _minkowski_sum_vrep_2d using ..HPolygonModule: HPolygon using LinearAlgebra: dot using Random: AbstractRNG, GLOBAL_RNG, shuffle @@ -13,7 +14,8 @@ using ReachabilityBase.Require: require @reexport import ..API: an_element, area, constraints_list, isoperationtype, rand, vertices_list, ∈, linear_map, permute, project, - σ, translate, translate! + σ, translate, translate!, convex_hull, intersection, + minkowski_sum @reexport import ..LazySets: remove_redundant_vertices, remove_redundant_vertices!, tohrep, tovrep import Base: convert @@ -35,6 +37,9 @@ include("permute.jl") include("project.jl") include("support_vector.jl") include("translate.jl") +include("convex_hull.jl") +include("intersection.jl") +include("minkowski_sum.jl") include("remove_redundant_vertices.jl") include("tohrep.jl") diff --git a/src/Sets/VPolygon/convex_hull.jl b/src/Sets/VPolygon/convex_hull.jl new file mode 100644 index 0000000000..433d4b0217 --- /dev/null +++ b/src/Sets/VPolygon/convex_hull.jl @@ -0,0 +1,26 @@ +""" + convex_hull(P::VPolygon, Q::VPolygon; [algorithm]::String="monotone_chain") + +Return the convex hull of two polygons in vertex representation. + +### Input + +- `P` -- polygon in vertex representation +- `Q` -- polygon in vertex representation +- `algorithm` -- (optional, default: "monotone_chain") the algorithm used to + compute the convex hull + +### Output + +A new polygon such that its vertices are the convex hull of the two polygons. + +### Notes + +The vertices of the output polygon are sorted in counter-clockwise fashion. +""" +function convex_hull(P::VPolygon, Q::VPolygon; + algorithm::String="monotone_chain") + vunion = [P.vertices; Q.vertices] + convex_hull!(vunion; algorithm=algorithm) + return VPolygon(vunion; apply_convex_hull=false) +end diff --git a/src/Sets/VPolygon/intersection.jl b/src/Sets/VPolygon/intersection.jl new file mode 100644 index 0000000000..5d2633e6a9 --- /dev/null +++ b/src/Sets/VPolygon/intersection.jl @@ -0,0 +1,37 @@ +""" + intersection(P1::VPolygon, P2::VPolygon; apply_convex_hull::Bool=true) + +Compute the intersection of two polygons in vertex representation. + +### Input + +- `P1` -- polygon in vertex representation +- `P2` -- polygon in vertex representation +- `apply_convex_hull` -- (default, optional: `true`) if `false`, skip the + computation of the convex hull of the resulting polygon + +### Output + +A `VPolygon`, or an `EmptySet` if the intersection is empty. + +### Algorithm + +This function applies the [Sutherland–Hodgman polygon clipping +algorithm](https://en.wikipedia.org/wiki/Sutherland%E2%80%93Hodgman_algorithm). +The implementation is based on the one found in +[rosetta code](http://www.rosettacode.org/wiki/Sutherland-Hodgman_polygon_clipping#Julia). +""" +function intersection(P1::VPolygon, P2::VPolygon; apply_convex_hull::Bool=true) + v1 = vertices_list(P1) + v2 = vertices_list(P2) + v12 = _intersection_vrep_2d(v1, v2) + + if isempty(v12) + require(@__MODULE__, :LazySets; fun_name="intersection") + + N = promote_type(eltype(P1), eltype(P2)) + return EmptySet{N}(2) + else + return VPolygon(v12; apply_convex_hull=apply_convex_hull) + end +end diff --git a/src/Sets/VPolygon/minkowski_sum.jl b/src/Sets/VPolygon/minkowski_sum.jl new file mode 100644 index 0000000000..3a2fc7eed7 --- /dev/null +++ b/src/Sets/VPolygon/minkowski_sum.jl @@ -0,0 +1,27 @@ +""" + minkowski_sum(P::VPolygon, Q::VPolygon) + +The Minkowski Sum of two polygons in vertex representation. + +### Input + +- `P` -- polygon in vertex representation +- `Q` -- polygon in vertex representation + +### Output + +A polygon in vertex representation. + +### Algorithm + +We treat each edge of the polygons as a vector, attaching them in polar order +(attaching the tail of the next vector to the head of the previous vector). The +resulting polygonal chain will be a polygon, which is the Minkowski sum of the +given polygons. This algorithm assumes that the vertices of `P` and `Q` are +sorted in counter-clockwise fashion and has linear complexity ``O(m+n)``, where +``m`` and ``n`` are the number of vertices of `P` and `Q`, respectively. +""" +function minkowski_sum(P::VPolygon, Q::VPolygon) + R = _minkowski_sum_vrep_2d(P.vertices, Q.vertices) + return VPolygon(R) +end diff --git a/src/Sets/VPolytope/VPolytopeModule.jl b/src/Sets/VPolytope/VPolytopeModule.jl index 77cdc1531f..3210702372 100644 --- a/src/Sets/VPolytope/VPolytopeModule.jl +++ b/src/Sets/VPolytope/VPolytopeModule.jl @@ -4,7 +4,8 @@ using Reexport, Requires using ..LazySets: AbstractPolytope, LazySet, LinearMapVRep, default_lp_solver, default_lp_solver_polyhedra, default_polyhedra_backend, - is_lp_infeasible, is_lp_optimal, linprog + is_lp_infeasible, is_lp_optimal, linprog, + _minkowski_sum_vrep_nd, _vertices_list using LinearAlgebra: dot using Random: AbstractRNG, GLOBAL_RNG using ReachabilityBase.Arrays: projection_matrix @@ -14,7 +15,8 @@ using ReachabilityBase.Require: require @reexport import ..API: constraints_list, dim, isoperationtype, rand, reflect, vertices_list, ∈, linear_map, permute, project, scale!, - ρ, σ, translate, translate! + ρ, σ, translate, translate!, cartesian_product, + convex_hull, minkowski_sum @reexport import ..LazySets: remove_redundant_vertices, tohrep, tovrep, _linear_map_vrep import Base: convert @@ -38,6 +40,9 @@ include("scale.jl") include("support_function.jl") include("support_vector.jl") include("translate.jl") +include("cartesian_product.jl") +include("convex_hull.jl") +include("minkowski_sum.jl") include("polyhedron.jl") include("remove_redundant_vertices.jl") diff --git a/src/Sets/VPolytope/cartesian_product.jl b/src/Sets/VPolytope/cartesian_product.jl new file mode 100644 index 0000000000..548016b8ac --- /dev/null +++ b/src/Sets/VPolytope/cartesian_product.jl @@ -0,0 +1,39 @@ +""" + cartesian_product(P1::VPolytope, P2::VPolytope; [backend]=nothing) + +Compute the Cartesian product of two polytopes in vertex representation. + +### Input + +- `P1` -- polytope in vertex representation +- `P2` -- polytope in vertex representation +- `backend` -- (optional, default: `nothing`) backend for polyhedral computation + +### Output + +The `VPolytope` obtained by the concrete Cartesian product of `P1` and `P2`. + +### Notes + +For further information on the supported backends see +[Polyhedra's documentation](https://juliapolyhedra.github.io/). +""" +function cartesian_product(P1::VPolytope, P2::VPolytope; backend=nothing) + require(@__MODULE__, :Polyhedra; fun_name="cartesian_product") + + return _cartesian_product_vrep(P1, P2; backend1=backend, backend2=backend) +end + +function _cartesian_product_vrep(P1, P2; backend1=nothing, backend2=nothing) + if isnothing(backend1) + backend1 = default_polyhedra_backend(P1) + end + if isnothing(backend2) + backend2 = default_polyhedra_backend(P2) + end + + P1′ = LazySets.polyhedron(P1; backend=backend1) + P2′ = LazySets.polyhedron(P2; backend=backend2) + Pout = Polyhedra.vcartesianproduct(P1′, P2′) + return VPolytope(Pout) +end diff --git a/src/Sets/VPolytope/convex_hull.jl b/src/Sets/VPolytope/convex_hull.jl new file mode 100644 index 0000000000..e0473d05cb --- /dev/null +++ b/src/Sets/VPolytope/convex_hull.jl @@ -0,0 +1,30 @@ +""" + convex_hull(P1::VPolytope, P2::VPolytope; [backend]=nothing) + +Compute the convex hull of two polytopes in vertex representation. + +### Input + +- `P1` -- polytope in vertex representation +- `P2` -- polytope in vertex representation +- `backend` -- (optional, default: `nothing`) the polyhedral computation backend + +### Output + +The `VPolytope` obtained by the concrete convex hull of `P1` and `P2`. + +### Notes + +This function takes the union of the vertices of each polytope and then relies +on a concrete convex-hull algorithm. +For low dimensions, a specialized implementation for polygons is used. +For higher dimensions, `convex_hull` relies on the polyhedral backend that can +be specified using the `backend` keyword argument. + +For performance reasons, it is suggested to use the `CDDLib.Library()` backend. +""" +function convex_hull(P1::VPolytope, P2::VPolytope; backend=nothing) + vunion = [P1.vertices; P2.vertices] + convex_hull!(vunion; backend=backend) + return VPolytope(vunion) +end diff --git a/src/Sets/VPolytope/minkowski_sum.jl b/src/Sets/VPolytope/minkowski_sum.jl new file mode 100644 index 0000000000..408283d413 --- /dev/null +++ b/src/Sets/VPolytope/minkowski_sum.jl @@ -0,0 +1,39 @@ +""" + minkowski_sum(P1::VPolytope, P2::VPolytope; + [apply_convex_hull]=true, + [backend]=nothing, + [solver]=nothing) + +Compute the Minkowski sum of two polytopes in vertex representation. + +### Input + +- `P1` -- polytope +- `P2` -- polytope +- `apply_convex_hull` -- (optional, default: `true`) if `true`, post-process the + pairwise sums using a convex-hull algorithm +- `backend` -- (optional, default: `nothing`) the backend for + polyhedral computations used to post-process with a + convex hull; see `default_polyhedra_backend(P1)` +- `solver` -- (optional, default: `nothing`) the backend used to + solve the linear program; see + `default_lp_solver_polyhedra(N)` + +### Output + +A new polytope in vertex representation whose vertices are the convex hull of +the sum of all possible sums of vertices of `P1` and `P2`. +""" +function minkowski_sum(P1::VPolytope, P2::VPolytope; + apply_convex_hull::Bool=true, backend=nothing, + solver=nothing) + @assert dim(P1) == dim(P2) "cannot compute the Minkowski sum of two " * + "polytopes of dimension $(dim(P1)) and $(dim(P2)), respectively" + + vlist1 = _vertices_list(P1, backend) + vlist2 = _vertices_list(P2, backend) + Vout = _minkowski_sum_vrep_nd(vlist1, vlist2; + apply_convex_hull=apply_convex_hull, + backend=backend, solver=solver) + return VPolytope(Vout) +end diff --git a/src/Sets/ZeroSet/ZeroSetModule.jl b/src/Sets/ZeroSet/ZeroSetModule.jl index 94e48f4794..823388904c 100644 --- a/src/Sets/ZeroSet/ZeroSetModule.jl +++ b/src/Sets/ZeroSet/ZeroSetModule.jl @@ -26,7 +26,6 @@ include("linear_map.jl") include("scale.jl") include("support_function.jl") include("translate.jl") - include("minkowski_sum.jl") include("element.jl") diff --git a/src/Utils/helper_functions.jl b/src/Utils/helper_functions.jl index 5800b86d54..774ada4871 100644 --- a/src/Utils/helper_functions.jl +++ b/src/Utils/helper_functions.jl @@ -269,3 +269,8 @@ function _is_linearcombination end # internal function; defined here due to dependency in submodules function _intersection_line2d end + +# internal functions; defined here due to dependency in submodules +function _minkowski_sum_hrep_preprocess end +function _minkowski_sum_vrep_2d end +function _minkowski_sum_vrep_nd end diff --git a/test/Sets/Polygon.jl b/test/Sets/Polygon.jl index f609279db1..35fbfd2767 100644 --- a/test/Sets/Polygon.jl +++ b/test/Sets/Polygon.jl @@ -230,6 +230,10 @@ for N in [Float64, Float32, Rational{Int}] oaux = VPolygon([N[0, 0], N[1 / 2, 0], N[0, 1 / 2]]) @test xaux ⊆ oaux && oaux ⊆ xaux # TODO use isequivalent @test LazySets._intersection_vrep_2d(paux.vertices, qaux.vertices) == xaux.vertices + # empty intersection + qaux = VPolygon([N[2, 2]]) + xaux = intersection(paux, qaux) + @test xaux == EmptySet{N}(2) # concrete intersection of H-rep and V-rep q1 = intersection(paux, p) diff --git a/test/Sets/Polytope.jl b/test/Sets/Polytope.jl index d4b2f7e1f8..36eda4d81d 100644 --- a/test/Sets/Polytope.jl +++ b/test/Sets/Polytope.jl @@ -578,10 +578,12 @@ for N in [Float64] @test !(BallInf(N[0, 0], N(1.01)) ⊆ P) # concrete minkowski sum - B = convert(VPolytope, BallInf(N[0, 0, 0], N(1))) - X = minkowski_sum(B, B) - twoB = 2.0 * B - @test X ⊆ twoB && twoB ⊆ X + for T in (HPolytope, VPolytope) + B = convert(T, BallInf(N[0, 0, 0], N(1))) + X = minkowski_sum(B, B) + twoB = 2.0 * B + @test X ⊆ twoB && twoB ⊆ X + end # concrete Reflection F4 = N[3 0; 0 7] diff --git a/test/Sets/SimpleSparsePolynomialZonotope.jl b/test/Sets/SimpleSparsePolynomialZonotope.jl index d22fceb38c..633a3785ff 100644 --- a/test/Sets/SimpleSparsePolynomialZonotope.jl +++ b/test/Sets/SimpleSparsePolynomialZonotope.jl @@ -70,6 +70,13 @@ for N in [Float64, Float32, Rational{Int}] @test genmat(CH1) == _g @test expmat(CH1) == _e + # convex hull with itself + CH2 = convex_hull(S, S) + Z1 = remove_redundant_generators(overapproximate(CH1, Zonotope)) + Z2 = remove_redundant_generators(overapproximate(CH2, Zonotope)) + @test_broken isequivalent(Z1, Z2) + @test Z1 ⊆ Z2 && center(Z1) == center(Z2) + S2 = SimpleSparsePolynomialZonotope(N[1 // 5, -3 // 5], N[1 0; 0 2//5], [1 0; 0 1]) Q = [N[0 1; 1 0], N[-16//5 6//5; 6//5 0]] diff --git a/test/Sets/SparsePolynomialZonotope.jl b/test/Sets/SparsePolynomialZonotope.jl index be7515dbdf..cf2145dcc5 100644 --- a/test/Sets/SparsePolynomialZonotope.jl +++ b/test/Sets/SparsePolynomialZonotope.jl @@ -43,6 +43,8 @@ for N in [Float64, Float32, Rational{Int}] @test genmat_indep(ESPZ) == hcat([1, 0]) @test expmat(ESPZ) == [1 0 3 1 0 1; 0 1 1 0 1 3] @test indexvector(ESPZ) == indexvector(PZ) + PZidx = SparsePolynomialZonotope(c, G, GI, E, 3:4) + @test_throws ArgumentError exact_sum(PZ, PZidx) TPZ = translate(PZ, N[1, 2]) @test center(TPZ) == N[5, 6]