Merge pull request #3368 from JuliaReach/schillic/1629

Add AbstractBallp interface

docs/src/lib/interfaces.md

@@ -553,3 +553,26 @@ dim(::AbstractPolynomialZonotope)
 * [Polynomial zonotope (DensePolynomialZonotope)](@ref def_DensePolynomialZonotope)
 * [Simplified sparse polynomial zonotope (SimpleSparsePolynomialZonotope)](@ref def_SimpleSparsePolynomialZonotope)
 * [Sparse polynomial zonotope (SparsePolynomialZonotope)](@ref def_SparsePolynomialZonotope)
+
+## [Balls in the p-norm (AbstractBallp)](@id def_AbstractBallp)
+
+A ball is a centrally-symmetric set with a characteristic p-norm.
+
+```@docs
+AbstractBallp
+```
+
+This interface defines the following functions:
+
+```@docs
+σ(::AbstractVector, ::AbstractBallp)
+ρ(::AbstractVector, ::AbstractBallp)
+∈(::AbstractVector, ::AbstractBallp)
+```
+
+### Implementations
+
+* [Ball1](@ref def_Ball1)
+* [Ball2](@ref def_Ball2)
+* [BallInf](@ref def_BallInf)
+* [Ballp](@ref def_Ballp)

docs/src/lib/sets/Ball1.md

@@ -6,11 +6,13 @@ CurrentModule = LazySets
 
 ```@docs
 Ball1
+center(::Ball1)
+radius_ball(::Ball1)
+ball_norm(::Ball1)
 σ(::AbstractVector, ::Ball1)
 ρ(::AbstractVector, ::Ball1)
 ∈(::AbstractVector, ::Ball1, ::Bool=false)
 vertices_list(::Ball1)
-center(::Ball1)
 rand(::Type{Ball1})
 constraints_list(::Ball1)
 translate(::Ball1, ::AbstractVector)

docs/src/lib/sets/Ball2.md

@@ -6,10 +6,12 @@ CurrentModule = LazySets
 
 ```@docs
 Ball2
+center(::Ball2)
+radius_ball(::Ball2)
+ball_norm(::Ball2)
 ρ(::AbstractVector, ::Ball2)
 σ(::AbstractVector, ::Ball2)
 ∈(::AbstractVector, ::Ball2)
-center(::Ball2)
 rand(::Type{Ball2})
 sample(::Ball2{N}, ::Int) where {N}
 translate(::Ball2, ::AbstractVector)

docs/src/lib/sets/BallInf.md

@@ -7,6 +7,8 @@ CurrentModule = LazySets
 ```@docs
 BallInf
 center(::BallInf)
+radius_ball(::BallInf)
+ball_norm(::BallInf)
 radius(::BallInf, ::Real=Inf)
 radius_hyperrectangle(::BallInf)
 radius_hyperrectangle(::BallInf, ::Int)

docs/src/lib/sets/Ballp.md

@@ -6,10 +6,9 @@ CurrentModule = LazySets
 
 ```@docs
 Ballp
-σ(::AbstractVector, ::Ballp)
-ρ(::AbstractVector, ::Ballp)
-∈(::AbstractVector, ::Ballp)
 center(::Ballp)
+radius_ball(::Ballp)
+ball_norm(::Ballp)
 rand(::Type{Ballp})
 translate(::Ballp, ::AbstractVector)
 translate!(::Ballp, ::AbstractVector)
@@ -31,3 +30,8 @@ Inherited from [`AbstractCentrallySymmetric`](@ref):
 * [`an_element`](@ref an_element(::AbstractCentrallySymmetric))
 * [`extrema`](@ref extrema(::AbstractCentrallySymmetric))
 * [`extrema`](@ref extrema(::AbstractCentrallySymmetric, ::Int))
+
+Inherited from [`AbstractBallp`](@ref):
+* [`σ`](@ref σ(::AbstractVector, ::AbstractBallp))
+* [`ρ`](@ref ρ(::AbstractVector, ::AbstractBallp))
+* [`∈`](@ref ∈(::AbstractVector, ::AbstractBallp))

src/ConcreteOperations/minkowski_sum.jl

@@ -507,3 +507,36 @@ function minkowski_sum(P1::SparsePolynomialZonotope,
     E = cat(expmat(P1), expmat(P2); dims=(1, 2))
     return SparsePolynomialZonotope(c, G, GI, E)
 end
+
+# Given two balls of the same p-norm, their Minkowski sum is again a p-norm ball,
+# which can be seen as follows:
+#
+# σ_Bp(d) = c + r \\frac{d}{‖d‖_p}
+#
+#   ρ_Bp(d)
+# = ⟨c + r \\frac{d}{‖d‖_p}, d⟩
+# = ⟨c, d⟩ + \\frac{r}{‖d‖_p} * ⟨d, d⟩
+#
+#   ρ_Bp¹⊕Bp²(d)
+# = ρ_Bp¹(d) + ρ_Bp²(d)
+# = ⟨c¹, d⟩ + \\frac{r¹}{‖d‖_p} * ⟨d, d⟩ + ⟨c², d⟩ + \\frac{r²}{‖d‖_p} * ⟨d, d⟩
+# = ⟨c¹ + c², d⟩ + \\frac{r¹ + r²}{‖d‖_p} * ⟨d, d⟩
+# = ρ_Bp³(d)
+#
+# where Bp³ = Ball(center(Bp¹) + center(Bp²), radius(Bp¹) + radius(Bp²))
+function minkowski_sum(B1::BT1, B2::BT2) where {BT<:AbstractBallp,BT1<:BT,BT2<:BT}
+    p = ball_norm(B1)
+    if ball_norm(B2) != p
+        throw(ArgumentError("this method only applies to balls of the same norm"))
+    end
+    return Ballp(p, center(B1) + center(B2), radius_ball(B1) + radius_ball(B2))
+end
+
+for B in (:Ball1, :BallInf)
+    @eval begin
+        # see AbstractBallp method
+        function minkowski_sum(B1::$B, B2::$B)
+            return $B(center(B1) + center(B2), radius_ball(B1) + radius_ball(B2))
+        end
+    end
+end

src/Interfaces/AbstractBallp.jl

@@ -0,0 +1,197 @@
+export AbstractBallp
+
+"""
+    AbstractBallp{N} <: AbstractCentrallySymmetric{N}
+
+Abstract type for p-norm balls.
+
+### Notes
+
+See [`Ballp`](@ref) for a standard implementation of this interface.
+
+Every concrete `AbstractBallp` must define the following methods:
+
+- `center(::AbstractBallp)` -- return the center
+- `radius_ball(::AbstractBallp)` -- return the ball radius
+- `ball_norm(::AbstractBallp)` -- return the characteristic norm
+
+The subtypes of `AbstractBallp`:
+
+```jldoctest; setup = :(using LazySets: subtypes)
+julia> subtypes(AbstractBallp)
+2-element Vector{Any}:
+ Ball2
+ Ballp
+```
+
+There are two further set types implementing the `AbstractBallp` interface, but
+they also implement other interfaces and hence cannot be subtypes: `Ball1` and
+`BallInf`.
+"""
+abstract type AbstractBallp{N} <: AbstractCentrallySymmetric{N} end
+
+function low(B::AbstractBallp)
+    return _low_AbstractBallp(B)
+end
+
+function _low_AbstractBallp(B::LazySet)
+    return center(B) .- radius_ball(B)
+end
+
+function low(B::AbstractBallp, i::Int)
+    return _low_AbstractBallp(B, i)
+end
+
+function _low_AbstractBallp(B::LazySet, i::Int)
+    return center(B, i) - radius_ball(B)
+end
+
+function high(B::AbstractBallp)
+    return _high_AbstractBallp(B)
+end
+
+function _high_AbstractBallp(B::LazySet)
+    return center(B) .+ radius_ball(B)
+end
+
+function high(B::AbstractBallp, i::Int)
+    return _high_AbstractBallp(B, i)
+end
+
+function _high_AbstractBallp(B::LazySet, i::Int)
+    return center(B, i) + radius_ball(B)
+end
+
+
+"""
+    σ(d::AbstractVector, B::AbstractBallp)
+
+Return the support vector of a ball in the p-norm in a given direction.
+
+### Input
+
+- `d` -- direction
+- `B` -- ball in the p-norm
+
+### Output
+
+The support vector in the given direction.
+If the direction has norm zero, the center of the ball is returned.
+
+### Algorithm
+
+The support vector of the unit ball in the ``p``-norm along direction ``d`` is:
+
+```math
+σ(d, \\mathcal{B}_p^n(0, 1)) = \\dfrac{\\tilde{v}}{‖\\tilde{v}‖_q},
+```
+where ``\\tilde{v}_i = \\frac{|d_i|^q}{d_i}`` if ``d_i ≠ 0`` and
+``\\tilde{v}_i = 0`` otherwise, for all ``i=1,…,n``, and ``q`` is the conjugate
+number of ``p``.
+By the affine transformation ``x = r\\tilde{x} + c``, one obtains that
+the support vector of ``\\mathcal{B}_p^n(c, r)`` is
+
+```math
+σ(d, \\mathcal{B}_p^n(c, r)) = \\dfrac{v}{‖v‖_q},
+```
+where ``v_i = c_i + r\\frac{|d_i|^q}{d_i}`` if ``d_i ≠ 0`` and ``v_i = 0``
+otherwise, for all ``i = 1, …, n``.
+"""
+function σ(d::AbstractVector, B::AbstractBallp)
+    p = ball_norm(B)
+    q = p / (p - 1)
+    v = similar(d)
+    N = promote_type(eltype(d), eltype(B))
+    @inbounds for (i, di) in enumerate(d)
+        v[i] = di == zero(N) ? di : abs.(di) .^ q / di
+    end
+    vnorm = norm(v, p)
+    if isapproxzero(vnorm)
+        svec = copy(center(B))
+    else
+        svec = center(B) .+ radius_ball(B) .* (v ./ vnorm)
+    end
+    return svec
+end
+
+"""
+    ρ(d::AbstractVector, B::AbstractBallp)
+
+Evaluate the support function of a ball in the p-norm in the given direction.
+
+### Input
+
+- `d` -- direction
+- `B` -- ball in the p-norm
+
+### Output
+
+Evaluation of the support function in the given direction.
+
+### Algorithm
+
+Let ``c`` and ``r`` be the center and radius of the ball ``B`` in the p-norm,
+respectively, and let ``q = \\frac{p}{p-1}``. Then:
+
+```math
+ρ(d, B) = ⟨d, c⟩ + r ‖d‖_q.
+```
+"""
+function ρ(d::AbstractVector, B::AbstractBallp)
+    p = ball_norm(B)
+    q = p / (p - 1)
+    return dot(d, center(B)) + radius_ball(B) * norm(d, q)
+end
+
+"""
+    ∈(x::AbstractVector, B::AbstractBallp)
+
+Check whether a given point is contained in a ball in the p-norm.
+
+### Input
+
+- `x` -- point/vector
+- `B` -- ball in the p-norm
+
+### Output
+
+`true` iff ``x ∈ B``.
+
+### Notes
+
+This implementation is worst-case optimized, i.e., it is optimistic and first
+computes (see below) the whole sum before comparing to the radius.
+In applications where the point is typically far away from the ball, a fail-fast
+implementation with interleaved comparisons could be more efficient.
+
+### Algorithm
+
+Let ``B`` be an ``n``-dimensional ball in the p-norm with radius ``r`` and let
+``c_i`` and ``x_i`` be the ball's center and the vector ``x`` in dimension
+``i``, respectively.
+Then ``x ∈ B`` iff ``\\left( ∑_{i=1}^n |c_i - x_i|^p \\right)^{1/p} ≤ r``.
+
+### Examples
+
+```jldoctest
+julia> B = Ballp(1.5, [1.0, 1.0], 1.)
+Ballp{Float64, Vector{Float64}}(1.5, [1.0, 1.0], 1.0)
+
+julia> [0.5, -0.5] ∈ B
+false
+
+julia> [0.5, 1.5] ∈ B
+true
+```
+"""
+function ∈(x::AbstractVector, B::AbstractBallp)
+    @assert length(x) == dim(B) "a vector of length $(length(x)) is " *
+                                "incompatible with a set of dimension $(dim(B))"
+    N = promote_type(eltype(x), eltype(B))
+    p = ball_norm(B)
+    sum = zero(N)
+    @inbounds for i in eachindex(x)
+        sum += abs(center(B, i) - x[i])^p
+    end
+    return _leq(sum^(one(N) / p), radius_ball(B))
+end

src/Interfaces/AbstractCentrallySymmetric.jl

@@ -21,9 +21,8 @@ The subtypes of `AbstractCentrallySymmetric`:
 
 ```jldoctest; setup = :(using LazySets: subtypes)
 julia> subtypes(AbstractCentrallySymmetric)
-3-element Vector{Any}:
- Ball2
- Ballp
+2-element Vector{Any}:
+ AbstractBallp
  Ellipsoid
 ```
 """

src/LazySets.jl

@@ -75,6 +75,7 @@ include("Interfaces/AbstractPolygon.jl")
 include("Interfaces/AbstractSingleton.jl")
 include("Interfaces/AbstractHPolygon.jl")
 include("Interfaces/AbstractAffineMap.jl")
+include("Interfaces/AbstractBallp.jl")
 
 # =============================
 # Types representing basic sets

src/Sets/Ball1.jl

@@ -74,6 +74,57 @@ function center(B::Ball1)
     return B.center
 end
 
+"""
+    radius_ball(B::Ball1)
+
+Return the ball radius of a ball in the 1-norm.
+
+### Input
+
+- `B` -- ball in the 1-norm
+
+### Output
+
+The ball radius.
+"""
+function radius_ball(B::Ball1)
+    return B.radius
+end
+
+"""
+    ball_norm(B::Ball1)
+
+Return the characteristic norm of a ball in the 1-norm.
+
+### Input
+
+- `B` -- ball in the 1-norm
+
+### Output
+
+The characteristic norm, which is `1`.
+"""
+function ball_norm(B::Ball1)
+    N = eltype(B)
+    return one(N)
+end
+
+function low(B::Ball1)
+    return _low_AbstractBallp(B)
+end
+
+function low(B::Ball1, i::Int)
+    return _low_AbstractBallp(B, i)
+end
+
+function high(B::Ball1)
+    return _high_AbstractBallp(B)
+end
+
+function high(B::Ball1, i::Int)
+    return _high_AbstractBallp(B, i)
+end
+
 """
     vertices_list(B::Ball1)