underapproximate with a Ball2 (#3055)

* fix kwarg name of chebyshev_center

* underapproximate with Ball2

* rename chebyshev_center to chebyshev_center_radius

docs/src/lib/interfaces.md

@@ -105,6 +105,7 @@ singleton_list(::ConvexSet)
 constraints(::ConvexSet)
 vertices(::ConvexSet)
 delaunay
+chebyshev_center_radius(::ConvexSet{N}) where {N}
 ```
 
 Plotting is available for general one- or two-dimensional `ConvexSet`s, provided
@@ -215,7 +216,6 @@ This interface defines the following functions:
 isuniversal(::AbstractPolyhedron{N}, ::Bool=false) where {N}
 constrained_dimensions(::AbstractPolyhedron)
 linear_map(::AbstractMatrix{NM}, ::AbstractPolyhedron{NP}) where {NM, NP}
-chebyshev_center(::AbstractPolyhedron{N}) where {N}
 an_element(::AbstractPolyhedron{N}) where {N}
 isbounded(::AbstractPolyhedron{N}) where {N}
 vertices_list(::AbstractPolyhedron)

docs/src/lib/sets/Ball2.md

@@ -14,7 +14,7 @@ rand(::Type{Ball2})
 sample(::Ball2{N}, ::Int) where {N}
 translate(::Ball2, ::AbstractVector)
 translate!(::Ball2, ::AbstractVector)
-chebyshev_center(::Ball2)
+chebyshev_center_radius(::Ball2)
 volume(::Ball2)
 ```
 Inherited from [`ConvexSet`](@ref):

docs/src/lib/sets/Interval.md

@@ -33,6 +33,8 @@ constraints_list(::Interval{N}) where {N}
 rectify(::Interval{N}) where {N}
 diameter(::Interval, ::Real=Inf)
 split(::Interval, ::AbstractVector{Int})
+chebyshev_center_radius(::Interval)
+fast_interval_pow(::IA.Interval, ::Int)
 ```
 Inherited from [`ConvexSet`](@ref):
 * [`singleton_list`](@ref singleton_list(::ConvexSet))

src/Approximations/underapproximate.jl

@@ -85,3 +85,26 @@ function underapproximate(X::ConvexSet{N}, ::Type{<:Hyperrectangle};
     solver = default_nln_solver(N)
     return _underapproximate_box(X, solver)
 end
+
+"""
+    underapproximate(X::LazySet, ::Type{<:Ball2})
+
+Compute the largest inscribed Euclidean ball in a set `X`.
+
+### Input
+
+- `X`     -- set
+- `Ball2` -- target type
+
+### Output
+
+A largest `Ball2` contained in `X`.
+
+### Algorithm
+
+We use `chebyshev_center_radius(X)`.
+"""
+function underapproximate(X::LazySet, ::Type{<:Ball2})
+    c, r = chebyshev_center_radius(X)
+    return Ball2(c, r)
+end

src/ConcreteOperations/difference.jl

@@ -3,8 +3,6 @@ export difference
 # alias for set difference
 import Base: \
 
-const IA = IntervalArithmetic
-
 """
     \\(X::ConvexSet, Y::ConvexSet)
 

src/Interfaces/AbstractPolyhedron_functions.jl

@@ -5,7 +5,6 @@ export constrained_dimensions,
        remove_redundant_constraints,
        remove_redundant_constraints!,
        linear_map,
-       chebyshev_center,
        an_element,
        vertices_list
 
@@ -852,56 +851,6 @@ function plot_recipe(P::AbstractPolyhedron{N}, ε=zero(N)) where {N}
     end
 end
 
-"""
-    chebyshev_center(P::AbstractPolyhedron{N};
-                     [get_radius]::Bool=false,
-                     [backend]=default_polyhedra_backend(P),
-                     [solver]=default_lp_solver_polyhedra(N; presolve=true)
-                     ) where {N}
-
-Compute the [Chebyshev center](https://en.wikipedia.org/wiki/Chebyshev_center)
-of a polytope.
-
-### Input
-
-- `P`          -- polytope
-- `get_radius` -- (optional; default: `false`) option to additionally return the
-                  radius of the largest ball enclosed by `P` around the
-                  Chebyshev center
-- `backend`    -- (optional; default: `default_polyhedra_backend(P)`) the
-                  backend for polyhedral computations
-- `solver`     -- (optional; default:
-                  `default_lp_solver_polyhedra(N; presolve=true)`) the LP
-                  solver passed to `Polyhedra`
-
-### Output
-
-If `get_radius` is `false`, the result is the Chebyshev center of `P`.
-If `get_radius` is `true`, the result is the pair `(c, r)` where `c` is the
-Chebyshev center of `P` and `r` is the radius of the largest ball with center
-`c` enclosed by `P`.
-
-### Notes
-
-The Chebyshev center is the center of a largest Euclidean ball enclosed by `P`.
-In general, the center of such a ball is not unique (but the radius is).
-"""
-function chebyshev_center(P::AbstractPolyhedron{N};
-                          get_radius::Bool=false,
-                          backend=default_polyhedra_backend(P),
-                          solver=default_lp_solver_polyhedra(N; presolve=true)
-                         ) where {N}
-    require(:Polyhedra; fun_name="chebyshev_center")
-    # convert to HPolyhedron to ensure `polyhedron` is applicable (see #1505)
-    Q = polyhedron(convert(HPolyhedron, P); backend=backend)
-    c, r = Polyhedra.chebyshevcenter(Q, solver)
-
-    if get_radius
-        return c, r
-    end
-    return c
-end
-
 """
     an_element(P::AbstractPolyhedron{N};
                [solver]=default_lp_solver(N)) where {N}

src/Interfaces/AbstractSingleton.jl

@@ -352,11 +352,8 @@ function ∈(S::AbstractSingleton, X::ConvexSet)
           "the implementations may differ)"))
 end
 
-function chebyshev_center(S::AbstractSingleton{N}; compute_radius::Bool=false) where {N}
-    if compute_radius
-        return element(S), zero(N)
-    end
-    return element(S)
+function chebyshev_center_radius(S::AbstractSingleton{N}) where {N}
+    return element(S), zero(N)
 end
 
 """

src/Interfaces/ConvexSet.jl

@@ -35,7 +35,8 @@ export ConvexSet,
        vertices,
        project,
        rectify,
-       permute
+       permute,
+       chebyshev_center_radius
 
 """
     ConvexSet{N} <: LazySet{N}
@@ -1488,3 +1489,46 @@ A new set corresponding to `X` where the dimensions have been permuted according
 to `p`.
 """
 function permute end
+
+"""
+    chebyshev_center_radius(P::ConvexSet{N};
+                            [backend]=default_polyhedra_backend(P),
+                            [solver]=default_lp_solver_polyhedra(N; presolve=true)
+                           ) where {N}
+
+Compute a [Chebyshev center](https://en.wikipedia.org/wiki/Chebyshev_center)
+and the corresponding radius of a polytopic set.
+
+### Input
+
+- `P`       -- polytopic set
+- `backend` -- (optional; default: `default_polyhedra_backend(P)`) the backend
+               for polyhedral computations
+- `solver`  -- (optional; default:
+               `default_lp_solver_polyhedra(N; presolve=true)`) the LP solver
+               passed to `Polyhedra`
+
+### Output
+
+The pair `(c, r)` where `c` is a Chebyshev center of `P` and `r` is the radius
+of the largest ball with center `c` enclosed by `P`.
+
+### Notes
+
+The Chebyshev center is the center of a largest Euclidean ball enclosed by `P`.
+In general, the center of such a ball is not unique, but the radius is.
+
+### Algorithm
+
+We call `Polyhedra.chebyshevcenter`.
+"""
+function chebyshev_center_radius(P::ConvexSet{N};
+                                 backend=default_polyhedra_backend(P),
+                                 solver=default_lp_solver_polyhedra(N; presolve=true)
+                                ) where {N}
+    require(:Polyhedra; fun_name="chebyshev_center")
+    # convert to HPolyhedron to ensure `polyhedron` is applicable (see #1505)
+    Q = polyhedron(convert(HPolyhedron, P); backend=backend)
+    c, r = Polyhedra.chebyshevcenter(Q, solver)
+    return c, r
+end

src/LazySets.jl

@@ -8,6 +8,7 @@ import GLPK, IntervalArithmetic, ReachabilityBase, JuMP, Random
 
 using IntervalArithmetic: AbstractInterval, mince
 import IntervalArithmetic: radius, ⊂
+const IA = IntervalArithmetic
 using LinearAlgebra: checksquare
 import LinearAlgebra: norm, ×, normalize, normalize!
 using Random: AbstractRNG, GLOBAL_RNG, SamplerType, shuffle, randperm

src/Sets/Ball2.jl

@@ -342,34 +342,27 @@ end
 
 
 """
-    chebyshev_center(B::Ball2; [compute_radius]::Bool=false)
+    chebyshev_center_radius(B::Ball2; [kwargs]...)
 
 Compute the [Chebyshev center](https://en.wikipedia.org/wiki/Chebyshev_center)
-of a ball in the 2-norm.
+and the corresponding radius of a ball in the 2-norm.
 
 ### Input
 
-- `B`              -- ball in the 2-norm
-- `compute_radius` -- (optional; default: `false`) option to additionally return
-                      the radius of the largest ball enclosed by `B` around the
-                      Chebyshev center
+- `B`      -- ball in the 2-norm
+- `kwargs` -- further keyword arguments (ignored)
 
 ### Output
 
-If `compute_radius` is `false`, the result is the Chebyshev center of `B`.
-If `compute_radius` is `true`, the result is the pair `(c, r)` where `c` is the
-Chebyshev center of `B` and `r` is the radius of the largest ball with center
-`c` enclosed by `B`.
+The pair `(c, r)` where `c` is the Chebyshev center of `B` and `r` is the radius
+of the largest ball with center `c` enclosed by `B`.
 
 ### Notes
 
 The Chebyshev center of a ball in the 2-norm is just the center of the ball.
 """
-function chebyshev_center(B::Ball2; compute_radius::Bool=false)
-    if compute_radius
-        return B.center, B.radius
-    end
-    return B.center
+function chebyshev_center_radius(B::Ball2; kwargs...)
+    return B.center, B.radius
 end
 
 """

src/Sets/Interval.jl

@@ -444,6 +444,10 @@ end
 # --- AbstractHyperrectangle interface functions ---
 
 
+function _radius(x::Interval{N}) where {N}
+    return (max(x) - min(x)) / N(2)
+end
+
 """
     radius_hyperrectangle(x::Interval{N}, i::Int) where {N}
 
@@ -460,7 +464,7 @@ The box radius in the given dimension.
 """
 function radius_hyperrectangle(x::Interval{N}, i::Int) where {N}
     @assert i == 1 "an interval is one-dimensional"
-    return (max(x) - min(x)) / N(2)
+    return _radius(x)
 end
 
 """
@@ -477,7 +481,7 @@ Return the box radius of an interval in every dimension.
 The box radius of the interval (a one-dimensional vector).
 """
 function radius_hyperrectangle(x::Interval)
-    return [radius_hyperrectangle(x, 1)]
+    return [_radius(x)]
 end
 
 """
@@ -702,11 +706,28 @@ function vertices_list(H::IntervalArithmetic.IntervalBox)
     return vertices_list(convert(Hyperrectangle, H))
 end
 
-function chebyshev_center(x::Interval{N}; compute_radius::Bool=false) where {N}
-    if compute_radius
-        return center(x), zero(N)
-    end
-    return center(x)
+"""
+    chebyshev_center_radius(x::Interval; [kwargs]...)
+
+Compute the [Chebyshev center](https://en.wikipedia.org/wiki/Chebyshev_center)
+and the corresponding radius of an interval.
+
+### Input
+
+- `x`      -- interval
+- `kwargs` -- further keyword arguments (ignored)
+
+### Output
+
+The pair `(c, r)` where `c` is the Chebyshev center of `x` and `r` is the radius
+of the largest ball with center `c` enclosed by `x`.
+
+### Notes
+
+The Chebyshev center of an interval is just the center of the interval.
+"""
+function chebyshev_center_radius(x::Interval; kwargs...)
+    return center(x), _radius(x)
 end
 
 """
@@ -730,7 +751,6 @@ use `a^n` from `IntervalArithmetic.jl`.
 
 Review after IntervalArithmetic.jl#388
 """
-const IA = IntervalArithmetic
 function fast_interval_pow(a::IA.Interval, n::Int)
     if iszero(n)
         return one(a)

test/Approximations/underapproximate.jl

@@ -7,3 +7,9 @@ for N in [Float64, Rational{Int}, Float32]
     U = underapproximate(X, Hyperrectangle)
     @test U ≈ Hyperrectangle(N[2, 2], ones(N, 2))
 end
+
+for N in [Float64, Float32]
+    X = VPolygon([N[1, 0], N[1, 2], N[-1, 2], N[-1, 1//3], N[-2//3, 0]])
+    U = underapproximate(X, Ball2)
+    @test U ≈ Ball2(N[0, 1], N(1))
+end

test/Sets/Ball2.jl

@@ -125,7 +125,8 @@ for N in [Float64, Float32]
     @test s isa AbstractVector{N} && s ∈ B
 
     # Chebyshev center
-    @test chebyshev_center(B) == center(B)
+    c, r = chebyshev_center_radius(B)
+    @test c == center(B) && r == B.radius
 
     # volume in dimension 2
     B = Ball2(zeros(N, 2), N(2))

test/Sets/Interval.jl

@@ -226,7 +226,6 @@ for N in [Float64, Float32, Rational{Int}]
     @test is_cyclic_permutation(vlistI, [SA[N(0)], SA[N(1)]])
 
     # Chebyshev center
-    c = chebyshev_center(x)
-    c2, r = chebyshev_center(x, compute_radius=true)
-    @test c == c2 == center(x) && r == zero(N)
+    c, r = chebyshev_center_radius(x)
+    @test c == center(x) && r == N(1//2)
 end

test/Sets/Polytope.jl

@@ -541,12 +541,11 @@ for N in [Float64]
 
         # Chebyshev center
         B = BallInf(N[0, 0], N(1))  # Chebyshev center is unique
-        c1 = chebyshev_center(B)
+        c1, r1 = chebyshev_center_radius(B)
         P = convert(HPolytope, B)
-        c2 = chebyshev_center(P)
-        c3, r = chebyshev_center(P; get_radius=true)
-        @test c1 == c2 == c3 == center(B) && c1 isa AbstractVector{N}
-        @test r == B.radius
+        c2, r2 = chebyshev_center_radius(P)
+        @test c1 == c2 == center(B) && c1 isa AbstractVector{N}
+        @test r1 == r2 == B.radius
 
         # concrete projection
         πP = project(P, [1])

test/Sets/Singleton.jl

@@ -156,7 +156,6 @@ for N in [Float64, Rational{Int}, Float32]
     @test permute(S, [2, 1, 3]) == Singleton(N[4, 3, 5])
 
     # Chebyshev center
-    c = chebyshev_center(S)
-    c2, r = chebyshev_center(S, compute_radius=true)
-    @test c == c2 == element(S) && r == zero(N)
+    c, r = chebyshev_center_radius(S)
+    @test c == element(S) && r == zero(N)
 end