diff --git a/docs/src/lib/interfaces.md b/docs/src/lib/interfaces.md index 8707db84e79..762624da912 100644 --- a/docs/src/lib/interfaces.md +++ b/docs/src/lib/interfaces.md @@ -562,6 +562,14 @@ A ball is a centrally-symmetric set with a characteristic p-norm. AbstractBallp ``` +This interface defines the following functions: + +```@docs +σ(::AbstractVector, ::AbstractBallp) +ρ(::AbstractVector, ::AbstractBallp) +∈(::AbstractVector, ::AbstractBallp) +``` + ### Implementations * [Ball1](@ref def_Ball1) diff --git a/docs/src/lib/sets/Ballp.md b/docs/src/lib/sets/Ballp.md index ec2b097060a..4ea6d218eab 100644 --- a/docs/src/lib/sets/Ballp.md +++ b/docs/src/lib/sets/Ballp.md @@ -9,9 +9,6 @@ Ballp center(::Ballp) radius_ball(::Ballp) ball_norm(::Ballp) -σ(::AbstractVector, ::Ballp) -ρ(::AbstractVector, ::Ballp) -∈(::AbstractVector, ::Ballp) rand(::Type{Ballp}) translate(::Ballp, ::AbstractVector) translate!(::Ballp, ::AbstractVector) @@ -33,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)) diff --git a/src/Interfaces/AbstractBallp.jl b/src/Interfaces/AbstractBallp.jl index 68ccaea804e..0b747976f03 100644 --- a/src/Interfaces/AbstractBallp.jl +++ b/src/Interfaces/AbstractBallp.jl @@ -61,3 +61,137 @@ 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 diff --git a/src/Sets/Ballp.jl b/src/Sets/Ballp.jl index b3fd40eca1c..1ad88a3a767 100644 --- a/src/Sets/Ballp.jl +++ b/src/Sets/Ballp.jl @@ -128,137 +128,6 @@ function ball_norm(B::Ballp) return B.p end -""" - σ(d::AbstractVector, B::Ballp) - -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::Ballp) - p = B.p - 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 = B.center - else - svec = @.(B.center + B.radius * (v / vnorm)) - end - return svec -end - -""" - ρ(d::AbstractVector, B::Ballp) - -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::Ballp) - q = B.p / (B.p - 1) - return dot(d, B.center) + B.radius * norm(d, q) -end - -""" - ∈(x::AbstractVector, B::Ballp) - -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::Ballp) - @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)) - sum = zero(N) - for i in eachindex(x) - sum += abs(B.center[i] - x[i])^B.p - end - return _leq(sum^(one(N) / B.p), B.radius) -end - """ rand(::Type{Ballp}; [N]::Type{<:Real}=Float64, [dim]::Int=2, [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing)