#3419 - Fast low/high/extrema for polyhedra

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

src/Interfaces/LazySet.jl

@@ -205,16 +205,38 @@ 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
 
 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
 
+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)
@@ -226,27 +248,73 @@ 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
 
 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
 
+function _high_vlist(X::LazySet)
+ n = dim(X)
+ vlist = vertices_list(X)
+ return [_high_vlist(vlist, i) for i in 1:n]
+end
+
 """
 ### Algorithm
 
 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)
 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
 
@@ -262,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, _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
 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
  # ===================================