[WIP] Evaluation distance to set

src/sets.jl

@@ -1,6 +1,7 @@
 # Sets
 
 # Note: When adding a new set, also add it to Utilities.Model.
+import LinearAlgebra
 
 """
     AbstractSet
@@ -80,6 +81,41 @@ function dual_set_type end
 
 dual_set_type(S::Type{<:AbstractSet}) = error("Dual type of $S is not implemented.")
 
+
+"""
+    AbstractDistance
+
+Distance function used to evaluate the distance from a point to a set.
+New subtypes of `AbstractDistance` must implement fallbacks for sets they don't cover and implement
+`distance_to_set(::Distance, v, s::S)` for sets they override the distance for.
+"""
+abstract type AbstractDistance end
+
+"""
+    DefaultDistance
+
+Default distance function, uses the Euclidean distance.
+"""
+struct DefaultDistance <: AbstractDistance end
+
+"""
+    distance_to_set(distance_definition, v, s)
+
+Compute the distance of a value to a set.
+When `v ∈ s`, the distance is zero (or all individual distances are zero).
+
+Each set `S` implements at least `distance_to_set(d::DefaultDistance, v::T, s::S)`
+with `T` of appropriate type for members of the set.
+"""
+function distance_to_set end
+
+distance_to_set(::AbstractDistance, v, s) = distance_to_set(DefaultDistance(), v, s)
+
+function _check_dimension(v::AbstractVector, s)
+    length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set"))
+    return nothing
+end
+
 """
     AbstractScalarSet
 
@@ -112,6 +148,11 @@ end
 dual_set(s::Reals) = Zeros(dimension(s))
 dual_set_type(::Type{Reals}) = Zeros
 
+function distance_to_set(::DefaultDistance, v::AbstractVector{T}, s::Reals) where {T <: Real}
+    _check_dimension(v, s)
+    return zero(T)
+end
+
 """
     Zeros(dimension)
 
@@ -124,6 +165,11 @@ end
 dual_set(s::Zeros) = Reals(dimension(s))
 dual_set_type(::Type{Zeros}) = Reals
 
+function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::Zeros)
+    _check_dimension(v, s)
+    return LinearAlgebra.norm2(v)
+end
+
 """
     Nonnegatives(dimension)
 
@@ -136,6 +182,11 @@ end
 dual_set(s::Nonnegatives) = copy(s)
 dual_set_type(::Type{Nonnegatives}) = Nonnegatives
 
+function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::Nonnegatives)
+    _check_dimension(v, s)
+    return LinearAlgebra.norm2(ifelse(vi < 0, -vi, zero(vi)) for vi in v)
+end
+
 """
     Nonpositives(dimension)
 
@@ -148,6 +199,11 @@ end
 dual_set(s::Nonpositives) = copy(s)
 dual_set_type(::Type{Nonpositives}) = Nonpositives
 
+function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::Nonpositives)
+    _check_dimension(v, s)
+    return LinearAlgebra.norm2(ifelse(vi > 0, vi, zero(vi)) for vi in v)
+end
+
 """
     GreaterThan{T <: Real}(lower::T)
 
@@ -179,6 +235,10 @@ function Base.:(==)(set1::S, set2::S) where S <: Union{GreaterThan, LessThan, Eq
     return constant(set1) == constant(set2)
 end
 
+distance_to_set(::DefaultDistance, v::Real, s::LessThan) = max(v - s.upper, zero(v))
+distance_to_set(::DefaultDistance, v::Real, s::GreaterThan) = max(s.lower - v, zero(v))
+distance_to_set(::DefaultDistance, v::Real, s::EqualTo) = abs(v - s.value)
+
 """
     Interval{T <: Real}(lower::T,upper::T)
 
@@ -206,6 +266,8 @@ Interval(s::LessThan{<:AbstractFloat}) = Interval(typemin(s.upper), s.upper)
 Interval(s::EqualTo{<:Real}) = Interval(s.value, s.value)
 Interval(s::Interval) = s
 
+distance_to_set(::DefaultDistance, v::T, s::Interval) where {T <: Real} = max(s.lower - v, v - s.upper, zero(T))
+
 """
     constant(s::Union{EqualTo, GreaterThan, LessThan})
 
@@ -227,6 +289,14 @@ end
 dual_set(s::NormInfinityCone) = NormOneCone(dimension(s))
 dual_set_type(::Type{NormInfinityCone}) = NormOneCone
 
+function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::NormInfinityCone)
+    _check_dimension(v, s)
+    t = v[1]
+    xs = v[2:end]
+    result = maximum(abs, xs) - t
+    return max(result, zero(result))
+end
+
 """
     NormOneCone(dimension)
 
@@ -239,6 +309,14 @@ end
 dual_set(s::NormOneCone) = NormInfinityCone(dimension(s))
 dual_set_type(::Type{NormOneCone}) = NormInfinityCone
 
+function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::NormOneCone)
+    _check_dimension(v, s)
+    t = v[1]
+    xs = v[2:end]
+    result = sum(abs, xs) - t
+    return max(result, zero(result))
+end
+
 """
     SecondOrderCone(dimension)
 
@@ -251,6 +329,14 @@ end
 dual_set(s::SecondOrderCone) = copy(s)
 dual_set_type(::Type{SecondOrderCone}) = SecondOrderCone
 
+function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::SecondOrderCone)
+    _check_dimension(v, s)
+    t = v[1]
+    xs = v[2:end]
+    result = LinearAlgebra.norm2(xs) - t
+    return max(result, zero(result))
+end
+
 """
     RotatedSecondOrderCone(dimension)
 
@@ -263,6 +349,16 @@ end
 dual_set(s::RotatedSecondOrderCone) = copy(s)
 dual_set_type(::Type{RotatedSecondOrderCone}) = RotatedSecondOrderCone
 
+function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::RotatedSecondOrderCone)
+    _check_dimension(v, s)
+    t = v[1]
+    u = v[2]
+    xs = v[3:end]
+    return LinearAlgebra.norm2(
+        (max(-t, zero(t)), max(-u, zero(u)), max(LinearAlgebra.dot(xs,xs) - 2 * t * u))
+    )
+end
+
 """
     GeometricMeanCone(dimension)
 
@@ -277,6 +373,22 @@ struct GeometricMeanCone <: AbstractVectorSet
     dimension::Int
 end
 
+function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::GeometricMeanCone)
+    _check_dimension(v, s)
+    t = v[1]
+    xs = v[2:end]
+    n = dimension(s) - 1
+    xresult = LinearAlgebra.norm2(
+        max.(-xs, zero(eltype(xs)))
+    )
+    # early returning if there exists x[i] < 0 to avoid complex sqrt
+    if xresult > 0
+        return xresult
+    end
+    result = t - prod(xs)^(inv(n))
+    return max(result, zero(result))
+end
+
 """
     ExponentialCone()
 
@@ -287,6 +399,17 @@ struct ExponentialCone <: AbstractVectorSet end
 dual_set(s::ExponentialCone) = DualExponentialCone()
 dual_set_type(::Type{ExponentialCone}) = DualExponentialCone
 
+function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::ExponentialCone)
+    _check_dimension(v, s)
+    x = v[1]
+    y = v[2]
+    z = v[3]
+    result = y * exp(x/y) - z
+    return LinearAlgebra.norm2(
+        (max(-y, zero(result)), max(result, zero(result)))
+    )
+end
+
 """
     DualExponentialCone()
 
@@ -297,6 +420,17 @@ struct DualExponentialCone <: AbstractVectorSet end
 dual_set(s::DualExponentialCone) = ExponentialCone()
 dual_set_type(::Type{DualExponentialCone}) = ExponentialCone
 
+function distance_to_set(::DefaultDistance, vs::AbstractVector{<:Real}, s::DualExponentialCone)
+    _check_dimension(vs, s)
+    u = vs[1]
+    v = vs[2]
+    w = vs[3]
+    result = -u*exp(v/u) - ℯ * w
+    return LinearAlgebra.norm2(
+        (max(u, zero(result)), max(result, zero(result)))
+    )
+end
+
 """
     PowerCone{T <: Real}(exponent::T)
 
@@ -309,6 +443,24 @@ end
 dual_set(s::PowerCone{T}) where T <: Real = DualPowerCone{T}(s.exponent)
 dual_set_type(::Type{PowerCone{T}}) where T <: Real = DualPowerCone{T}
 
+function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::PowerCone)
+    _check_dimension(v, s)
+    x = v[1]
+    y = v[2]
+    z = v[3]
+    e = s.exponent
+    # early return to avoid complex exponent results
+    if x < 0 || y < 0
+       return LinearAlgebra.norm2(
+            (max(-x, zero(x)), max(-y, zero(x)))
+        )
+    end
+    result = abs(z) - x^e * y^(1-e)
+    return LinearAlgebra.norm2(
+        max(result, zero(result))
+    )
+end
+
 """
     DualPowerCone{T <: Real}(exponent::T)
 
@@ -321,6 +473,19 @@ end
 dual_set(s::DualPowerCone{T}) where T <: Real = PowerCone{T}(s.exponent)
 dual_set_type(::Type{DualPowerCone{T}}) where T <: Real = PowerCone{T}
 
+function distance_to_set(::DefaultDistance, vs::AbstractVector{<:Real}, s::DualPowerCone)
+    _check_dimension(vs, s)
+    u = vs[1]
+    v = vs[2]
+    w = vs[3]
+    e = s.exponent
+    ce = 1-e
+    result = abs(w) - (u/e)^e * (v/ce)^ce
+    return LinearAlgebra.norm2(
+        (max(-u, zero(result)), max(-v, zero(result)), max(result, zero(result)))
+    )
+end
+
 dimension(s::Union{ExponentialCone, DualExponentialCone, PowerCone, DualPowerCone}) = 3
 
 function Base.:(==)(set1::S, set2::S) where S <: Union{PowerCone, DualPowerCone}
@@ -342,6 +507,20 @@ struct RelativeEntropyCone <: AbstractVectorSet
     dimension::Int
 end
 
+function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, set::RelativeEntropyCone)
+    _check_dimension(v, s)
+    n = (dimension(set)-1) ÷ 2
+    u = v[1]
+    v = v[2:(n+1)]
+    w = v[(n+2):end]
+    s = sum(w[i] * log(w[i]/v[i]) for i in eachindex(w))
+    result = s - u
+    return LinearAlgebra.norm2(push!(
+        max.(v[2:end], zero(result)),
+        max(result, zero(result)),
+    ))
+end
+
 """
     NormSpectralCone(row_dim, column_dim)
 
@@ -356,6 +535,15 @@ end
 dual_set(s::NormSpectralCone) = NormNuclearCone(s.row_dim, s.column_dim)
 dual_set_type(::Type{NormSpectralCone}) = NormNuclearCone
 
+function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::NormSpectralCone)
+    _check_dimension(v, s)
+    t = v[1]
+    m = reshape(v[2:end], (s.row_dim, s.column_dim))
+    s1 = LinearAlgebra.svd(m).S[1]
+    result = s1 - t
+    return max(result, zero(result))
+end
+
 """
     NormNuclearCone(row_dim, column_dim)
 
@@ -370,6 +558,15 @@ end
 dual_set(s::NormNuclearCone) = NormSpectralCone(s.row_dim, s.column_dim)
 dual_set_type(::Type{NormNuclearCone}) = NormSpectralCone
 
+function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::NormNuclearCone)
+    _check_dimension(v, s)
+    t = v[1]
+    m = reshape(v[2:end], (s.row_dim, s.column_dim))
+    s1 = sum(LinearAlgebra.svd(m).S)
+    result = s1 - t
+    return max(result, zero(result))
+end
+
 dimension(s::Union{NormSpectralCone, NormNuclearCone}) = 1 + s.row_dim * s.column_dim
 
 """
@@ -653,6 +850,14 @@ The set ``\\{ 0, 1 \\}``.
 """
 struct ZeroOne <: AbstractScalarSet end
 
+function distance_to_set(::DefaultDistance, v::T, ::ZeroOne) where {T <: Real}
+    return min(abs(v - zero(T)), abs(v - one(T)))
+end
+
+function distance_to_set(::DefaultDistance, v::Real, ::Integer)
+    return min(abs(v - floor(v)), abs(v - ceil(v)))
+end
+
 """
     Semicontinuous{T <: Real}(lower::T,upper::T)
 
@@ -690,6 +895,13 @@ struct SOS1{T <: Real} <: AbstractVectorSet
     weights::Vector{T}
 end
 
+# return the element-wise distance to zero, with the greatest element to 0
+function distance_to_set(::DefaultDistance, v::AbstractVector{T}, ::SOS1) where {T <: Real}
+    _check_dimension(v, s)
+    _, i = findmax(abs.(v))
+    return LinearAlgebra.norm2([v[j] for j = eachindex(v) if j != i])
+end
+
 """
     SOS2{T <: Real}(weights::Vector{T})
 
@@ -755,6 +967,19 @@ end
 
 dimension(::IndicatorSet) = 2
 
+# takes in input [z, f(x)]
+function distance_to_set(d::DefaultDistance, v::AbstractVector{T}, s::IndicatorSet{A}) where {A, T <: Real}
+    _check_dimension(v, s)
+    z = v[1]
+    # inactive constraint
+    if A === ACTIVATE_ON_ONE && isapprox(z, 0) || A === ACTIVATE_ON_ZERO && isapprox(z, 1)
+        return zeros(T, 2)
+    end
+    return LinearAlgebra.norm2(
+        (distance_to_set(d, z, ZeroOne()), distance_to_set(v[2], s.set))
+    )
+end
+
 function Base.copy(set::IndicatorSet{A,S}) where {A,S}
     return IndicatorSet{A}(copy(set.set))
 end
@@ -852,6 +1077,7 @@ gives the `n-1`-dimensional nonnegative orthant. However
 function supports_dimension_update(::Type{<:AbstractVectorSet})
     return false
 end
+
 function supports_dimension_update(::Type{<:Union{
     Reals, Zeros, Nonnegatives, Nonpositives}})
     return true