Skip to content

Commit

Permalink
generalize code from Ballp to AbstractBallp
Browse files Browse the repository at this point in the history
  • Loading branch information
schillic committed Aug 3, 2023
1 parent 726487f commit 6d8e7f3
Show file tree
Hide file tree
Showing 4 changed files with 147 additions and 134 deletions.
8 changes: 8 additions & 0 deletions docs/src/lib/interfaces.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
8 changes: 5 additions & 3 deletions docs/src/lib/sets/Ballp.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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))
134 changes: 134 additions & 0 deletions src/Interfaces/AbstractBallp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
131 changes: 0 additions & 131 deletions src/Sets/Ballp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down

0 comments on commit 6d8e7f3

Please sign in to comment.