Merge pull request #3617 from JuliaReach/schillic/3419

#3419 - Fast `low`/`high`/`extrema` for polyhedra

src/Interfaces/AbstractPolyhedron_functions.jl

@@ -1248,3 +1248,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

@@ -204,16 +204,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)
@@ -225,27 +247,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
 
@@ -261,6 +329,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
     # ===================================