better low/high/extrema for VPolygon/VPolytope

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}
 

src/Sets/VPolygon/VPolygonModule.jl

@@ -3,19 +3,20 @@ 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
+ is_right_turn, _area_vlist, _extrema_vlist, _high_vlist,
+ _intersection_vrep_2d, _linear_map_vrep, _low_vlist,
+ _minkowski_sum_vrep_2d
 using ..HPolygonModule: HPolygon
 using LinearAlgebra: dot
 using Random: AbstractRNG, GLOBAL_RNG, shuffle
 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")

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

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

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

src/Sets/VPolytope/VPolytopeModule.jl

@@ -4,19 +4,19 @@ 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
 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")

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

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

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

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
  # ===================================

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))