[WIP] Evaluation distance to set

src/MathOptInterface.jl

@@ -113,7 +113,6 @@ include("error.jl")
 include("indextypes.jl")
 include("functions.jl")
 include("sets.jl")
-include("attributes.jl")
 include("constraints.jl")
 include("modifications.jl")
 include("variables.jl")

src/sets.jl

@@ -80,6 +80,22 @@ function dual_set_type end
 
 dual_set_type(S::Type{<:AbstractSet}) = error("Dual type of $S is not implemented.")
 
+"""
+    distance_to_set(v, s)
+
+Compute the distance of a value to a set.
+For some vector-valued sets, can return a vector of distances.
+When `v ∈ s`, the distance is zero (or all individual distances are zero).
+
+Each set `S` implements `distance_to_set(v::T, s::S)` with `T` of appropriate type.
+"""
+function distance_to_set end
+
+function _check_dimension(v::AbstractVector, s)
+    length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set"))
+    return nothing
+end
+
 """
     AbstractScalarSet
 
@@ -112,6 +128,11 @@ end
 dual_set(s::Reals) = Zeros(dimension(s))
 dual_set_type(::Type{Reals}) = Zeros
 
+function distance_to_set(v::AbstractVector{<:Real}, s::Reals)
+    _check_dimension(v, s)
+    return zeros(length(v))
+end
+
 """
     Zeros(dimension)
 
@@ -124,6 +145,11 @@ end
 dual_set(s::Zeros) = Reals(dimension(s))
 dual_set_type(::Type{Zeros}) = Reals
 
+function distance_to_set(v::AbstractVector{<:Real}, s::Zeros)
+    _check_dimension(v, s)
+    return abs.(v)
+end
+
 """
     Nonnegatives(dimension)
 
@@ -136,6 +162,11 @@ end
 dual_set(s::Nonnegatives) = copy(s)
 dual_set_type(::Type{Nonnegatives}) = Nonnegatives
 
+function distance_to_set(v::AbstractVector{<:Real}, s::Nonnegatives)
+    _check_dimension(v, s)
+    return [ifelse(vi < 0, -vi, zero(vi)) for vi in v]
+end
+
 """
     Nonpositives(dimension)
 
@@ -148,6 +179,11 @@ end
 dual_set(s::Nonpositives) = copy(s)
 dual_set_type(::Type{Nonpositives}) = Nonpositives
 
+function distance_to_set(v::AbstractVector{<:Real}, s::Nonpositives)
+    _check_dimension(v, s)
+    return [ifelse(vi > 0, vi, zero(vi)) for vi in v]
+end
+
 """
     GreaterThan{T <: Real}(lower::T)
 
@@ -179,6 +215,10 @@ function Base.:(==)(set1::S, set2::S) where S <: Union{GreaterThan, LessThan, Eq
     return constant(set1) == constant(set2)
 end
 
+distance_to_set(v::Real, s::LessThan) = max(v - s.upper, zero(v))
+distance_to_set(v::Real, s::GreaterThan) = max(s.lower - v, zero(v))
+distance_to_set(v::Real, s::EqualTo) = abs(v - s.value)
+
 """
     Interval{T <: Real}(lower::T,upper::T)
 
@@ -206,6 +246,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(v::T, s::Interval) where {T <: Real} = max(s.lower - v, v - s.upper, zero(T))
+
 """
     constant(s::Union{EqualTo, GreaterThan, LessThan})
 
@@ -239,6 +281,14 @@ end
 dual_set(s::NormOneCone) = NormInfinityCone(dimension(s))
 dual_set_type(::Type{NormOneCone}) = NormInfinityCone
 
+function distance_to_set(v::AbstractVector{<:Real}, s::NormOneCone)
+    _check_dimension(v, s)
+    t = v[1]
+    xs = v[2:end]
+    result = t - sum(abs, xs)
+    return max(result, zero(result))
+end
+
 """
     SecondOrderCone(dimension)
 
@@ -251,6 +301,14 @@ end
 dual_set(s::SecondOrderCone) = copy(s)
 dual_set_type(::Type{SecondOrderCone}) = SecondOrderCone
 
+function distance_to_set(v::AbstractVector{<:Real}, s::SecondOrderCone)
+    _check_dimension(v, s)
+    t = v[1]
+    xs = v[2:end]
+    result = LinearAlgebra.norm(xs) - t
+    return max(result, zero(result)) # avoids sqrt
+end
+
 """
     RotatedSecondOrderCone(dimension)
 
@@ -263,6 +321,18 @@ end
 dual_set(s::RotatedSecondOrderCone) = copy(s)
 dual_set_type(::Type{RotatedSecondOrderCone}) = RotatedSecondOrderCone
 
+function distance_to_set(v::AbstractVector{<:Real}, s::RotatedSecondOrderCone)
+    _check_dimension(v, s)
+    t = v[1]
+    u = v[2]
+    xs = v[3:end]
+    return [
+        max(-t, zero(t)),
+        max(-u, zero(u)),
+        max(dot(xs,xs) - 2 * t * u),
+    ]
+end
+
 """
     GeometricMeanCone(dimension)
 
@@ -272,6 +342,18 @@ struct GeometricMeanCone <: AbstractVectorSet
     dimension::Int
 end
 
+function distance_to_set(v::AbstractVector{<:Real}, s::GeometricMeanCone)
+    _check_dimension(v, s)
+    t = v[1]
+    xs = v[2:end]
+    n = dimension(s) - 1
+    result = t^n - prod(xs)
+    return push!(
+        max.(xs, zero(result)), # x >= 0
+        max(result, zero(result)),
+    )
+end
+
 """
     ExponentialCone()
 
@@ -282,6 +364,15 @@ struct ExponentialCone <: AbstractVectorSet end
 dual_set(s::ExponentialCone) = DualExponentialCone()
 dual_set_type(::Type{ExponentialCone}) = DualExponentialCone
 
+function distance_to_set(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 [max(y, zero(result)), max(result, zero(result))]
+end
+
 """
     DualExponentialCone()
 
@@ -292,6 +383,15 @@ struct DualExponentialCone <: AbstractVectorSet end
 dual_set(s::DualExponentialCone) = ExponentialCone()
 dual_set_type(::Type{DualExponentialCone}) = ExponentialCone
 
+function distance_to_set(v::AbstractVector{<:Real}, s::DualExponentialCone)
+    _check_dimension(v, s)
+    u = v[1]
+    v = v[2]
+    w = v[3]
+    result = -u*exp(v/u) - ℯ * w
+    return [max(-u, zero(result)), max(result, zero(result))]
+end
+
 """
     PowerCone{T <: Real}(exponent::T)
 
@@ -304,6 +404,16 @@ 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(v::AbstractVector{<:Real}, s::PowerCone)
+    _check_dimension(v, s)
+    x = v[1]
+    y = v[2]
+    z = v[3]
+    e = s.exponent
+    result = abs(z) - x^e * y^(1-e)
+    return [max(-x, zero(result)), max(-y, zero(result)), max(result, zero(result))]
+end
+
 """
     DualPowerCone{T <: Real}(exponent::T)
 
@@ -316,6 +426,17 @@ 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(v::AbstractVector{<:Real}, s::DualPowerCone)
+    _check_dimension(v, s)
+    u = v[1]
+    v = v[2]
+    w = v[3]
+    e = s.exponent
+    ce = 1-e
+    result = abs(w) - (u/e)^e * (v/ce)^ce
+    return [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}
@@ -331,6 +452,20 @@ struct RelativeEntropyCone <: AbstractVectorSet
     dimension::Int
 end
 
+function distance_to_set(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 push!(
+        max.(v[2:end], zero(result)),
+        max(result, zero(result)),
+    )
+end
+
 """
     NormSpectralCone(row_dim, column_dim)
 
@@ -345,6 +480,15 @@ end
 dual_set(s::NormSpectralCone) = NormNuclearCone(s.row_dim, s.column_dim)
 dual_set_type(::Type{NormSpectralCone}) = NormNuclearCone
 
+function distance_to_set(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)
 
@@ -359,6 +503,15 @@ end
 dual_set(s::NormNuclearCone) = NormSpectralCone(s.row_dim, s.column_dim)
 dual_set_type(::Type{NormNuclearCone}) = NormSpectralCone
 
+function distance_to_set(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
 
 """
@@ -642,6 +795,14 @@ The set ``\\{ 0, 1 \\}``.
 """
 struct ZeroOne <: AbstractScalarSet end
 
+function distance_to_set(v::T, ::ZeroOne) where {T <: Real}
+    return min(abs(v - zero(T)), abs(v - one(T)))
+end
+
+function distance_to_set(v::Real, ::Integer)
+    return min(abs(v - floor(v)), abs(v - ceil(v)))
+end
+
 """
     Semicontinuous{T <: Real}(lower::T,upper::T)
 
@@ -679,6 +840,23 @@ 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(v::AbstractVector{T}, ::SOS1) where {T <: Real}
+    _check_dimension(v, s)
+    d = distance_to_set(v, Zeros(length(v)))
+    m = maximum(d)
+    if m ≈ zero(T)
+        return d
+    end
+    # removing greatest distance
+    for i in eachindex(d)
+        @inbounds if d[i] == m
+            d[i] = zero(T)
+            return d
+        end
+    end
+end
+
 """
     SOS2{T <: Real}(weights::Vector{T})
 
@@ -744,6 +922,17 @@ end
 
 dimension(::IndicatorSet) = 2
 
+# takes in input [z, f(x)]
+function distance_to_set(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 [distance_to_set(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