From a7e84e8a9c0e75d610cd339c6649241f51681c8e Mon Sep 17 00:00:00 2001 From: Christian Schilling Date: Mon, 22 Aug 2022 11:47:00 +0200 Subject: [PATCH] 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 | 2 +- docs/src/lib/sets/Ball2.md | 2 +- docs/src/lib/sets/Interval.md | 2 + src/Approximations/underapproximate.jl | 23 +++++++++ src/ConcreteOperations/difference.jl | 2 - .../AbstractPolyhedron_functions.jl | 51 ------------------- src/Interfaces/AbstractSingleton.jl | 7 +-- src/Interfaces/ConvexSet.jl | 46 ++++++++++++++++- src/LazySets.jl | 1 + src/Sets/Ball2.jl | 23 +++------ src/Sets/Interval.jl | 36 ++++++++++--- test/Approximations/underapproximate.jl | 6 +++ test/Sets/Ball2.jl | 3 +- test/Sets/Interval.jl | 5 +- test/Sets/Polytope.jl | 9 ++-- test/Sets/Singleton.jl | 5 +- 16 files changed, 127 insertions(+), 96 deletions(-) diff --git a/docs/src/lib/interfaces.md b/docs/src/lib/interfaces.md index ba6c636dde..dbebb0975b 100644 --- a/docs/src/lib/interfaces.md +++ b/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) diff --git a/docs/src/lib/sets/Ball2.md b/docs/src/lib/sets/Ball2.md index d7d9fa8b8a..c00528808e 100644 --- a/docs/src/lib/sets/Ball2.md +++ b/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): diff --git a/docs/src/lib/sets/Interval.md b/docs/src/lib/sets/Interval.md index 00ea3b0c11..89ea99a334 100644 --- a/docs/src/lib/sets/Interval.md +++ b/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)) diff --git a/src/Approximations/underapproximate.jl b/src/Approximations/underapproximate.jl index 60012ddad0..cece9e973c 100644 --- a/src/Approximations/underapproximate.jl +++ b/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 diff --git a/src/ConcreteOperations/difference.jl b/src/ConcreteOperations/difference.jl index 7fa3a49bfb..15bf21e2c6 100644 --- a/src/ConcreteOperations/difference.jl +++ b/src/ConcreteOperations/difference.jl @@ -3,8 +3,6 @@ export difference # alias for set difference import Base: \ -const IA = IntervalArithmetic - """ \\(X::ConvexSet, Y::ConvexSet) diff --git a/src/Interfaces/AbstractPolyhedron_functions.jl b/src/Interfaces/AbstractPolyhedron_functions.jl index 94f84a966c..45ab051df1 100644 --- a/src/Interfaces/AbstractPolyhedron_functions.jl +++ b/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} diff --git a/src/Interfaces/AbstractSingleton.jl b/src/Interfaces/AbstractSingleton.jl index 27988a1ccc..23aa77db7e 100644 --- a/src/Interfaces/AbstractSingleton.jl +++ b/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 """ diff --git a/src/Interfaces/ConvexSet.jl b/src/Interfaces/ConvexSet.jl index c524160d83..7e93429383 100644 --- a/src/Interfaces/ConvexSet.jl +++ b/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 diff --git a/src/LazySets.jl b/src/LazySets.jl index 01c370367e..ce4f6350b7 100644 --- a/src/LazySets.jl +++ b/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 diff --git a/src/Sets/Ball2.jl b/src/Sets/Ball2.jl index dce66e9980..b59aa5e2c2 100644 --- a/src/Sets/Ball2.jl +++ b/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 """ diff --git a/src/Sets/Interval.jl b/src/Sets/Interval.jl index d923f36cb4..bfb87dfd82 100644 --- a/src/Sets/Interval.jl +++ b/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) diff --git a/test/Approximations/underapproximate.jl b/test/Approximations/underapproximate.jl index 4f3274148c..0530da3b3a 100644 --- a/test/Approximations/underapproximate.jl +++ b/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 diff --git a/test/Sets/Ball2.jl b/test/Sets/Ball2.jl index c8b6958199..d4100d5789 100644 --- a/test/Sets/Ball2.jl +++ b/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)) diff --git a/test/Sets/Interval.jl b/test/Sets/Interval.jl index 46d1051463..b391ceb53d 100644 --- a/test/Sets/Interval.jl +++ b/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 diff --git a/test/Sets/Polytope.jl b/test/Sets/Polytope.jl index f0ca562929..76fb54e982 100644 --- a/test/Sets/Polytope.jl +++ b/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]) diff --git a/test/Sets/Singleton.jl b/test/Sets/Singleton.jl index 03169bfed9..cb66f449df 100644 --- a/test/Sets/Singleton.jl +++ b/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