From df13c741eceecf8c2c67d55ad07222fd2de4b305 Mon Sep 17 00:00:00 2001 From: schillic Date: Fri, 26 Jul 2024 09:41:44 +0200 Subject: [PATCH 1/2] fast low/high/extrema for polyhedra --- .../AbstractPolyhedron_functions.jl | 91 +++++++++++++++++++ src/Interfaces/LazySet.jl | 12 +++ test/Sets/Polyhedron.jl | 12 +++ 3 files changed, 115 insertions(+) diff --git a/src/Interfaces/AbstractPolyhedron_functions.jl b/src/Interfaces/AbstractPolyhedron_functions.jl index ec04e9c785..aa3d5ecc70 100644 --- a/src/Interfaces/AbstractPolyhedron_functions.jl +++ b/src/Interfaces/AbstractPolyhedron_functions.jl @@ -1249,3 +1249,94 @@ Determine whether a polyhedron is equivalent to a hyperplane. ``a·x ≤ b`` and ``-a·x ≤ -b``. """ function ishyperplanar(::AbstractPolyhedron) end + +function extrema(P::AbstractPolyhedron) + if dim(P) == 1 + l, h = _extrema_1d(P) + return ([l], [h]) + end + return _extrema_lowhigh(P) +end + +function extrema(P::AbstractPolyhedron, i::Int) + if dim(P) == 1 + @assert i == 1 "invalid index $i for a set of 1 dimension" + return _extrema_1d(P) + end + return _extrema_lowhigh(P, i) +end + +function _extrema_1d(P::AbstractPolyhedron{N}) where {N} + l = N(-Inf) + h = N(Inf) + @inbounds for c in constraints(P) + a = c.a[1] + if a > zero(N) + # a·x <= b => x <= b/a + h = min(h, c.b / a) + else + # HalfSpace `0·x <= b` is not allowed (type constraint) + @assert a < zero(N) "invalid constraint" + + # -a·x <= -b => x >= b/a + l = max(l, c.b / a) + end + end + return (l, h) +end + +function low(P::AbstractPolyhedron) + if dim(P) == 1 + l = _low_1d(P) + return [l] + end + return _low(P) +end + +function low(P::AbstractPolyhedron, i::Int) + if dim(P) == 1 + @assert i == 1 "invalid index $i for a set of 1 dimension" + return _low_1d(P) + end + return _low(P, i) +end + +function _low_1d(P::AbstractPolyhedron{N}) where {N} + l = N(-Inf) + @inbounds for c in constraints(P) + a = c.a[1] + if a < zero(N) + # -a·x <= -b => x >= b/a + l = max(l, c.b / a) + end + end + return l +end + +function high(P::AbstractPolyhedron) + if dim(P) == 1 + l = _high_1d(P) + return [l] + end + return _high(P) +end + +function high(P::AbstractPolyhedron, i::Int) + if dim(P) == 1 + @assert i == 1 "invalid index $i for a set of 1 dimension" + return _high_1d(P) + end + return _high(P, i) +end + +function _high_1d(P::AbstractPolyhedron{N}) where {N} + h = N(Inf) + @inbounds for c in constraints(P) + a = c.a[1] + if a > zero(N) + # a·x <= b => x <= b/a + h = min(h, c.b / a) + end + end + return h +end diff --git a/src/Interfaces/LazySet.jl b/src/Interfaces/LazySet.jl index a3a49cf70d..703071e831 100644 --- a/src/Interfaces/LazySet.jl +++ b/src/Interfaces/LazySet.jl @@ -211,6 +211,10 @@ end The default implementation applies `low` in each dimension. """ function low(X::LazySet) + return _low(X) +end + +function _low(X::LazySet) n = dim(X) return [low(X, i) for i in 1:n] end @@ -232,6 +236,10 @@ end The default implementation applies `high` in each dimension. """ function high(X::LazySet) + return _high(X) +end + +function _high(X::LazySet) n = dim(X) return [high(X, i) for i in 1:n] end @@ -242,6 +250,10 @@ end The default implementation computes the extrema via `low` and `high`. """ function extrema(X::LazySet, i::Int) + return _extrema_lowhigh(X, i) +end + +function _extrema_lowhigh(X::LazySet, i::Int) l = low(X, i) h = high(X, i) return (l, h) diff --git a/test/Sets/Polyhedron.jl b/test/Sets/Polyhedron.jl index 0ebe2996ce..bf21d8997b 100644 --- a/test/Sets/Polyhedron.jl +++ b/test/Sets/Polyhedron.jl @@ -173,6 +173,18 @@ for N in [Float64, Rational{Int}, Float32] @test permute(P, 1:2) == P @test permute(P, [2, 1]) == HPolyhedron(HalfSpace[HalfSpace(N[2, 1], N(3)), HalfSpace(N[-2, -1], N(-3))]) + + # low/high/extrema + P = HPolyhedron(HalfSpace[HalfSpace(N[1, 0], N(3)), HalfSpace(N[-1, 0], N(-2))]) + @test low(P) == N[2, -Inf] && low(P, 1) == N(2) && low(P, 2) == N(-Inf) + @test high(P) == N[3, Inf] && high(P, 1) == N(3) && high(P, 2) == N(Inf) + @test extrema(P) == (N[2, -Inf], N[3, Inf]) && extrema(P, 1) == (N(2), N(3)) && + extrema(P, 2) == (N(-Inf), N(Inf)) + # 1D + P = HPolyhedron(HalfSpace[HalfSpace(N[2], N(6)), HalfSpace(N[-1], N(1))]) + @test low(P) == N[-1] && low(P, 1) == N(-1) + @test high(P) == N[3] && high(P, 1) == N(3) + @test extrema(P) == (N[-1], N[3]) && extrema(P, 1) == (N(-1), N(3)) end # default Float64 constructors From db761156f1c177b41d127b274414e3122b8bfa08 Mon Sep 17 00:00:00 2001 From: schillic Date: Fri, 26 Jul 2024 11:05:13 +0200 Subject: [PATCH 2/2] better low/high/extrema for VPolygon/VPolytope --- src/Interfaces/LazySet.jl | 80 +++++++++++++++++++++++++++ src/Sets/VPolygon/VPolygonModule.jl | 16 ++++-- src/Sets/VPolygon/extrema.jl | 8 +++ src/Sets/VPolygon/high.jl | 8 +++ src/Sets/VPolygon/low.jl | 8 +++ src/Sets/VPolytope/VPolytopeModule.jl | 17 +++--- src/Sets/VPolytope/extrema.jl | 8 +++ src/Sets/VPolytope/high.jl | 8 +++ src/Sets/VPolytope/low.jl | 8 +++ test/Sets/Polygon.jl | 19 +++++++ test/Sets/Polytope.jl | 21 +++++++ 11 files changed, 188 insertions(+), 13 deletions(-) create mode 100644 src/Sets/VPolygon/extrema.jl create mode 100644 src/Sets/VPolygon/high.jl create mode 100644 src/Sets/VPolygon/low.jl create mode 100644 src/Sets/VPolytope/extrema.jl create mode 100644 src/Sets/VPolytope/high.jl create mode 100644 src/Sets/VPolytope/low.jl diff --git a/src/Interfaces/LazySet.jl b/src/Interfaces/LazySet.jl index 703071e831..a729c2f98c 100644 --- a/src/Interfaces/LazySet.jl +++ b/src/Interfaces/LazySet.jl @@ -205,6 +205,18 @@ function _low(X::LazySet{N}, i::Int) where {N} return -ρ(d, X) end +function _low_vlist(X::LazySet, i::Int) + return _low_vlist(vertices_list(X), i) +end + +function _low_vlist(vlist::AbstractVector{<:AbstractVector{N}}, i::Int) where {N} + l = N(Inf) + @inbounds for v in vlist + l = min(l, v[i]) + end + return l +end + """ ### Algorithm @@ -219,6 +231,12 @@ function _low(X::LazySet) return [low(X, i) for i in 1:n] end +function _low_vlist(X::LazySet) + n = dim(X) + vlist = vertices_list(X) + return [_low_vlist(vlist, i) for i in 1:n] +end + # Note: this method cannot be documented due to a bug in Julia function high(X::LazySet, i::Int) return _high(X, i) @@ -230,6 +248,18 @@ function _high(X::LazySet{N}, i::Int) where {N} return ρ(d, X) end +function _high_vlist(X::LazySet, i::Int) + return _high_vlist(vertices_list(X), i) +end + +function _high_vlist(vlist::AbstractVector{<:AbstractVector{N}}, i::Int) where {N} + h = N(-Inf) + @inbounds for v in vlist + h = max(h, v[i]) + end + return h +end + """ ### Algorithm @@ -244,6 +274,12 @@ function _high(X::LazySet) return [high(X, i) for i in 1:n] end +function _high_vlist(X::LazySet) + n = dim(X) + vlist = vertices_list(X) + return [_high_vlist(vlist, i) for i in 1:n] +end + """ ### Algorithm @@ -259,6 +295,26 @@ function _extrema_lowhigh(X::LazySet, i::Int) return (l, h) end +function _extrema_vlist(X::LazySet, i::Int) + return _extrema_vlist(vertices_list(X), i) +end + +function _extrema_vlist(vlist::AbstractVector{<:AbstractVector{N}}, i::Int) where {N} + if isempty(vlist) + return (N(Inf), N(-Inf)) + end + l = h = @inbounds vlist[1][i] + @inbounds for v in vlist + vi = v[i] + if vi > h + h = vi + elseif vi < l + l = vi + end + end + return (l, h) +end + """ ### Algorithm @@ -274,6 +330,30 @@ function _extrema_lowhigh(X::LazySet) return (l, h) end +function _extrema_vlist(X::LazySet{N}) where {N} + n = dim(X) + if n <= 0 + return (N[], N[]) + end + vlist = vertices_list(X) + if isempty(vlist) + return (fill(N(Inf), n), fill(N(-Inf), n)) + end + l = copy(@inbounds vlist[1]) + h = copy(l) + @inbounds for v in vlist + for i in 1:n + vi = v[i] + if vi > h[i] + h[i] = vi + elseif vi < l[i] + l[i] = vi + end + end + end + return (l, h) +end + """ plot_vlist(X::S, ε::Real) where {S<:LazySet} diff --git a/src/Sets/VPolygon/VPolygonModule.jl b/src/Sets/VPolygon/VPolygonModule.jl index 6292f5e708..d31f48f3d7 100644 --- a/src/Sets/VPolygon/VPolygonModule.jl +++ b/src/Sets/VPolygon/VPolygonModule.jl @@ -3,8 +3,9 @@ module VPolygonModule using Reexport, Requires using ..LazySets: AbstractPolygon, LazySet, AbstractHPolygon, halfspace_left, - is_right_turn, _area_vlist, _intersection_vrep_2d, - _linear_map_vrep, _minkowski_sum_vrep_2d, _to_colVector + is_right_turn, _area_vlist, _extrema_vlist, _high_vlist, + _intersection_vrep_2d, _linear_map_vrep, _low_vlist, + _minkowski_sum_vrep_2d, _to_colVector using ..HPolygonModule: HPolygon using LinearAlgebra: dot using Random: AbstractRNG, GLOBAL_RNG, shuffle @@ -12,10 +13,10 @@ using ReachabilityBase.Arrays: isabove, rand_pos_neg_zerosum_vector using ReachabilityBase.Distribution: reseed! using ReachabilityBase.Require: require -@reexport import ..API: an_element, area, constraints_list, isoperationtype, - rand, vertices_list, ∈, linear_map, permute, project, - σ, translate, translate!, convex_hull, intersection, - minkowski_sum +@reexport import ..API: an_element, area, constraints_list, extrema, high, + isoperationtype, low, rand, vertices_list, ∈, + linear_map, permute, project, σ, translate, translate!, + convex_hull, intersection, minkowski_sum @reexport import ..LazySets: remove_redundant_vertices, remove_redundant_vertices!, tohrep, tovrep import Base: convert @@ -28,7 +29,10 @@ include("VPolygon.jl") include("an_element.jl") include("area.jl") include("constraints_list.jl") +include("extrema.jl") +include("high.jl") include("isoperationtype.jl") +include("low.jl") include("rand.jl") include("vertices_list.jl") include("in.jl") diff --git a/src/Sets/VPolygon/extrema.jl b/src/Sets/VPolygon/extrema.jl new file mode 100644 index 0000000000..b555e5dc69 --- /dev/null +++ b/src/Sets/VPolygon/extrema.jl @@ -0,0 +1,8 @@ +function extrema(P::VPolygon) + return _extrema_vlist(P) +end + +function extrema(P::VPolygon, i::Int) + @assert 1 <= i <= dim(P) "invalid index $i for set of dimension $(dim(P))" + return _extrema_vlist(P, i) +end diff --git a/src/Sets/VPolygon/high.jl b/src/Sets/VPolygon/high.jl new file mode 100644 index 0000000000..0db798912e --- /dev/null +++ b/src/Sets/VPolygon/high.jl @@ -0,0 +1,8 @@ +function high(P::VPolygon) + return _high_vlist(P) +end + +function high(P::VPolygon, i::Int) + @assert 1 <= i <= dim(P) "invalid index $i for set of dimension $(dim(P))" + return _high_vlist(P, i) +end diff --git a/src/Sets/VPolygon/low.jl b/src/Sets/VPolygon/low.jl new file mode 100644 index 0000000000..f9f4f6442f --- /dev/null +++ b/src/Sets/VPolygon/low.jl @@ -0,0 +1,8 @@ +function low(P::VPolygon) + return _low_vlist(P) +end + +function low(P::VPolygon, i::Int) + @assert 1 <= i <= dim(P) "invalid index $i for set of dimension $(dim(P))" + return _low_vlist(P, i) +end diff --git a/src/Sets/VPolytope/VPolytopeModule.jl b/src/Sets/VPolytope/VPolytopeModule.jl index 019231b5b2..cc907bc7a5 100644 --- a/src/Sets/VPolytope/VPolytopeModule.jl +++ b/src/Sets/VPolytope/VPolytopeModule.jl @@ -4,8 +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, - _minkowski_sum_vrep_nd + is_lp_infeasible, is_lp_optimal, linprog, _extrema_vlist, + _high_vlist, _low_vlist, _minkowski_sum_vrep_nd using LinearAlgebra: dot using Random: AbstractRNG, GLOBAL_RNG using ReachabilityBase.Arrays: projection_matrix @@ -13,10 +13,10 @@ using ReachabilityBase.Comparison: _ztol using ReachabilityBase.Distribution: reseed! using ReachabilityBase.Require: require -@reexport import ..API: constraints_list, dim, isoperationtype, rand, reflect, - vertices_list, ∈, linear_map, permute, project, scale!, - ρ, σ, translate, translate!, cartesian_product, - convex_hull, minkowski_sum +@reexport import ..API: constraints_list, dim, extrema, high, isoperationtype, + low, rand, reflect, vertices_list, ∈, linear_map, + permute, project, scale!, ρ, σ, translate, translate!, + cartesian_product, convex_hull, minkowski_sum @reexport import ..LazySets: remove_redundant_vertices, tohrep, tovrep import ..LazySets: _linear_map_vrep import Base: convert @@ -28,11 +28,14 @@ include("VPolytope.jl") include("constraints_list.jl") include("dim.jl") +include("extrema.jl") +include("high.jl") +include("isoperationtype.jl") +include("low.jl") include("rand.jl") include("reflect.jl") include("vertices_list.jl") include("in.jl") -include("isoperationtype.jl") include("linear_map.jl") include("permute.jl") include("project.jl") diff --git a/src/Sets/VPolytope/extrema.jl b/src/Sets/VPolytope/extrema.jl new file mode 100644 index 0000000000..adaac28913 --- /dev/null +++ b/src/Sets/VPolytope/extrema.jl @@ -0,0 +1,8 @@ +function extrema(P::VPolytope) + return _extrema_vlist(P) +end + +function extrema(P::VPolytope, i::Int) + @assert 1 <= i <= dim(P) "invalid index $i for set of dimension $(dim(P))" + return _extrema_vlist(P, i) +end diff --git a/src/Sets/VPolytope/high.jl b/src/Sets/VPolytope/high.jl new file mode 100644 index 0000000000..1d5779e465 --- /dev/null +++ b/src/Sets/VPolytope/high.jl @@ -0,0 +1,8 @@ +function high(P::VPolytope) + return _high_vlist(P) +end + +function high(P::VPolytope, i::Int) + @assert 1 <= i <= dim(P) "invalid index $i for set of dimension $(dim(P))" + return _high_vlist(P, i) +end diff --git a/src/Sets/VPolytope/low.jl b/src/Sets/VPolytope/low.jl new file mode 100644 index 0000000000..9bd79296aa --- /dev/null +++ b/src/Sets/VPolytope/low.jl @@ -0,0 +1,8 @@ +function low(P::VPolytope) + return _low_vlist(P) +end + +function low(P::VPolytope, i::Int) + @assert 1 <= i <= dim(P) "invalid index $i for set of dimension $(dim(P))" + return _low_vlist(P, i) +end diff --git a/test/Sets/Polygon.jl b/test/Sets/Polygon.jl index 2ccf07e36e..f64c492a6f 100644 --- a/test/Sets/Polygon.jl +++ b/test/Sets/Polygon.jl @@ -206,6 +206,25 @@ for N in [Float64, Float32, Rational{Int}] N[9, 12], N[7, 12], N[4, 10]]) end + # low/high/extrema + p2 = VPolygon([N[5, 1], N[4, 0], N[3, 1], N[4, 2]]) + @test low(p2) == N[3, 0] && low(p2, 1) == N(3) && low(p2, 2) == N(0) + @test high(p2) == N[5, 2] && high(p2, 1) == N(5) && high(p2, 2) == N(2) + @test extrema(p2) == (N[3, 0], N[5, 2]) && extrema(p2, 1) == (N(3), N(5)) && + extrema(p2, 2) == (N(0), N(2)) + # singleton + p2 = VPolygon([N[1, 2]]) + @test low(p2) == N[1, 2] && low(p2, 1) == N(1) && low(p2, 2) == N(2) + @test high(p2) == N[1, 2] && high(p2, 1) == N(1) && high(p2, 2) == N(2) + @test extrema(p2) == (N[1, 2], N[1, 2]) && extrema(p2, 1) == (N(1), N(1)) && + extrema(p2, 2) == (N(2), N(2)) + # empty polygon + p2 = VPolygon{N}() + @test low(p2) == N[Inf, Inf] && low(p2, 1) == N(Inf) && low(p2, 2) == N(Inf) + @test high(p2) == N[-Inf, -Inf] && high(p2, 1) == N(-Inf) && high(p2, 2) == N(-Inf) + @test extrema(p2) == (N[Inf, Inf], N[-Inf, -Inf]) && extrema(p2, 1) == (N(Inf), N(-Inf)) && + extrema(p2, 2) == (N(Inf), N(-Inf)) + # =================================== # Concrete intersection # =================================== diff --git a/test/Sets/Polytope.jl b/test/Sets/Polytope.jl index ccc82122b5..646de47125 100644 --- a/test/Sets/Polytope.jl +++ b/test/Sets/Polytope.jl @@ -230,6 +230,27 @@ for N in [Float64, Rational{Int}, Float32] # vertices_list function @test vertices_list(p) == p.vertices + # low/high/extrema + p2 = VPolytope([N[5, 1], N[4, 0], N[3, 1], N[4, 2]]) + @test low(p2) == N[3, 0] && low(p2, 1) == N(3) && low(p2, 2) == N(0) + @test high(p2) == N[5, 2] && high(p2, 1) == N(5) && high(p2, 2) == N(2) + @test extrema(p2) == (N[3, 0], N[5, 2]) && extrema(p2, 1) == (N(3), N(5)) && + extrema(p2, 2) == (N(0), N(2)) + # singleton + p2 = VPolytope([N[1, 2]]) + @test low(p2) == N[1, 2] && low(p2, 1) == N(1) && low(p2, 2) == N(2) + @test high(p2) == N[1, 2] && high(p2, 1) == N(1) && high(p2, 2) == N(2) + @test extrema(p2) == (N[1, 2], N[1, 2]) && extrema(p2, 1) == (N(1), N(1)) && + extrema(p2, 2) == (N(2), N(2)) + # empty polytope cannot determine dimension and returns empty vectors + p2 = VPolytope{N}() + @test low(p2) == N[] + @test_throws AssertionError low(p2, 1) + @test high(p2) == N[] + @test_throws AssertionError high(p2, 1) + @test extrema(p2) == (N[], N[]) + @test_throws AssertionError extrema(p2, 1) + if test_suite_polyhedra V = VPolytope(polyhedron(p))