[Utilities] add more distance_to_set implementations

src/Utilities/distance_to_set.jl

@@ -213,3 +213,174 @@ function distance_to_set(
     # Projection to the point (t, x) + 0.5 * (|x|_2 - t, (t/|x|_2 - 1) * x)
     return sqrt(2) / 2 * abs(t - rhs)
 end
+
+function distance_to_set(
+    ::ProjectionUpperBoundDistance,
+    x::AbstractVector{T},
+    set::MOI.RotatedSecondOrderCone,
+) where {T<:Real}
+    _check_dimension(x, set)
+    t, u, xs = x[1], x[2], @view(x[3:end])
+    element_distance = (
+        min(t, zero(T)),  # t >= 0
+        min(u, zero(T)),  # u >= 0
+        max(LinearAlgebra.dot(xs, xs) - 2 * t * u, zero(T)),
+    )
+    return LinearAlgebra.norm(element_distance, 2)
+end
+
+function distance_to_set(
+    ::ProjectionUpperBoundDistance,
+    x::AbstractVector{T},
+    set::MOI.ExponentialCone,
+) where {T<:Real}
+    _check_dimension(x, set)
+    if x[2] <= 0  # Project to x[2] = 1
+        element_distance = (x[2] - one(T), max(exp(x[1]) - x[3], zero(T)))
+        return LinearAlgebra.norm(element_distance, 2)
+    end
+    return max(x[2] * exp(x[1] / x[2]) - x[3], zero(T))
+end
+
+function distance_to_set(
+    ::ProjectionUpperBoundDistance,
+    x::AbstractVector{T},
+    set::MOI.DualExponentialCone,
+) where {T<:Real}
+    _check_dimension(x, set)
+    if x[1] >= 0  # Project to x[1] = -1
+        element_distance =
+            (x[1] - -one(T), max(exp(-x[2]) - exp(1) * x[3], zero(T)))
+        return LinearAlgebra.norm(element_distance, 2)
+    end
+    return max(-x[1] * exp(x[2] / x[1]) - exp(1) * x[3], zero(T))
+end
+
+function distance_to_set(
+    ::ProjectionUpperBoundDistance,
+    x::AbstractVector{T},
+    set::MOI.GeometricMeanCone,
+) where {T<:Real}
+    _check_dimension(x, set)
+    t, xs = x[1], @view(x[2:end])
+    if any(<(zero(T)), xs)  # Project to x = 0
+        return LinearAlgebra.norm((min.(xs, zero(T)), max(t, zero(T))), 2)
+    end
+    return max(t - prod(xs)^inv(MOI.dimension(set) - 1), zero(T))
+end
+
+function distance_to_set(
+    ::ProjectionUpperBoundDistance,
+    x::AbstractVector{T},
+    set::MOI.PowerCone,
+) where {T<:Real}
+    _check_dimension(x, set)
+    α = set.exponent
+    if x[1] < 0 || x[2] < 0  # Project to x = 0
+        return LinearAlgebra.norm(
+            (min(x[1], zero(T)), min(x[2], zero(T)), x[3]),
+            2,
+        )
+    end
+    return max(abs(x[3]) - x[1]^α * x[2]^(1 - α), zero(T))
+end
+
+function distance_to_set(
+    ::ProjectionUpperBoundDistance,
+    x::AbstractVector{T},
+    set::MOI.DualPowerCone,
+) where {T<:Real}
+    _check_dimension(x, set)
+    α = set.exponent
+    if x[1] < 0 || x[2] < 0  # Project to x = 0
+        return LinearAlgebra.norm(
+            (min(x[1], zero(T)), min(x[2], zero(T)), x[3]),
+            2,
+        )
+    end
+    return max(abs(x[3]) - (x[1] / α)^α * (x[2] / (1 - α))^(1 - α), zero(T))
+end
+
+function distance_to_set(
+    ::ProjectionUpperBoundDistance,
+    x::AbstractVector{T},
+    set::MOI.NormOneCone,
+) where {T<:Real}
+    _check_dimension(x, set)
+    return max(sum(abs, @view(x[2:end])) - x[1], zero(T))
+end
+
+function distance_to_set(
+    ::ProjectionUpperBoundDistance,
+    x::AbstractVector{T},
+    set::MOI.NormInfinityCone,
+) where {T<:Real}
+    _check_dimension(x, set)
+    return max(maximum(abs, @view(x[2:end])) - x[1], zero(T))
+end
+
+function distance_to_set(
+    ::ProjectionUpperBoundDistance,
+    x::AbstractVector{T},
+    set::MOI.RelativeEntropyCone,
+) where {T<:Real}
+    _check_dimension(x, set)
+    n = div(MOI.dimension(set) - 1, 2)
+    u, v, w = x[1], @view(x[2:(n+1)]), @view(x[(n+2):end])
+    _to_one(x) = x <= zero(T) ? one(T) : x
+    if any(<=(zero(T)), v) || any(<=(zero(T)), w)  # Project to v = w = 1
+        v_p, w_p = _to_one.(v), _to_one.(w)
+        element_distance = (
+            v_p .- v,
+            w_p .- w,
+            max(
+                sum(w_p[i] * log(w_p[i] / v_p[i]) for i in eachindex(w)) - u,
+                zero(T),
+            ),
+        )
+        return LinearAlgebra.norm(element_distance, 2)
+    end
+    return max(sum(w[i] * log(w[i] / v[i]) for i in eachindex(w)) - u, zero(T))
+end
+
+function distance_to_set(
+    ::ProjectionUpperBoundDistance,
+    x::AbstractVector{T},
+    set::MOI.HyperRectangle,
+) where {T<:Real}
+    _check_dimension(x, set)
+    element_distance =
+        (max.(set.lower .- x, zero(T)), max.(x .- set.upper, zero(T)))
+    return LinearAlgebra.norm(element_distance, 2)
+end
+
+function distance_to_set(
+    ::ProjectionUpperBoundDistance,
+    x::AbstractVector{T},
+    set::MOI.NormCone,
+) where {T<:Real}
+    _check_dimension(x, set)
+    return max(LinearAlgebra.norm(@view(x[2:end]), set.p) - x[1], zero(T))
+end
+
+function distance_to_set(
+    ::ProjectionUpperBoundDistance,
+    x::AbstractVector{T},
+    set::MOI.SOS1,
+) where {T<:Real}
+    _check_dimension(x, set)
+    _, i = findmax(abs, x)
+    return LinearAlgebra.norm2([x[j] for j in eachindex(x) if j != i])
+end
+
+function distance_to_set(
+    ::ProjectionUpperBoundDistance,
+    x::AbstractVector{T},
+    set::MOI.SOS2,
+) where {T<:Real}
+    _check_dimension(x, set)
+    p = sortperm(set.weights)
+    pairs = collect(zip(p[1:end-1], p[2:end]))
+    _, k = findmax(abs(x[i]) + abs(x[j]) for (i, j) in pairs)
+    return LinearAlgebra.norm2([x[i] for i in eachindex(x) if !(i in pairs[k])])
+end

test/Utilities/distance_to_set.jl

@@ -22,6 +22,19 @@ function runtests()
     return
 end
 
+function _test_set(set, pairs...; mismatch = nothing)
+    if mismatch !== nothing
+        @test_throws(
+            DimensionMismatch,
+            MOI.Utilities.distance_to_set(mismatch, set),
+        )
+    end
+    for (x, d) in pairs
+        @test MOI.Utilities.distance_to_set(x, set) ≈ d
+    end
+    return
+end
+
 function test_unsupported()
     @test_throws(
         ErrorException,
@@ -141,6 +154,165 @@ function test_secondordercone()
     return
 end
 
+function test_rotatedsecondordercone()
+    _test_set(
+        MOI.RotatedSecondOrderCone(4),
+        [1.0, 1.0, 1.0, 1.0] => 0.0,
+        [-1.0, 1.0, 1.0, 1.0] => sqrt(1 + 4^2);
+        mismatch = [1.0],
+    )
+    return
+end
+
+function test_exponential()
+    _test_set(
+        MOI.ExponentialCone(),
+        [1.0, 1.0, 1.0] => exp(1) - 1,
+        [1.0, 1.0, 3.0] => 0.0,
+        [1.0, -1.0, 1.0] => sqrt(2^2 + (exp(1) - 1)^2);
+        mismatch = [1.0],
+    )
+    return
+end
+
+function test_dualexponential()
+    _test_set(
+        MOI.DualExponentialCone(),
+        [1.0, 1.0, 1.0] => 2.0,
+        [-1.0, 1.0, 3.0] => 0.0,
+        [-2.0, 3.0, 0.1] => 2 * exp(3 / -2) - 0.1 * exp(1);
+        mismatch = [1.0],
+    )
+    return
+end
+
+function test_geometricmeancone()
+    _test_set(
+        MOI.GeometricMeanCone(3),
+        [1.0, 1.0, 1.0] => 0.0,
+        [1.5, 1.0, 2.0] => 1.5 - sqrt(2),
+        [3.5, 3.0, 2.0] => 3.5 - sqrt(6),
+        [1.5, -1.0, 2.0] => sqrt(1 + 1.5^2);
+        mismatch = [1.0],
+    )
+    return
+end
+
+function test_powercone()
+    _test_set(
+        MOI.PowerCone(0.5),
+        [1.0, 1.0, 1.0] => 0.0,
+        [-1.0, 1.0, 1.0] => sqrt(2),
+        [1.0, -1.0, 2.0] => sqrt(5),
+        [1.5, 1.0, -2.0] => 2 - 1.5^0.5 * 1^0.5,
+        [1.5, 1.0, 2.0] => 2 - 1.5^0.5 * 1^0.5;
+        mismatch = [1.0],
+    )
+    return
+end
+
+function test_dualpowercone()
+    _test_set(
+        MOI.DualPowerCone(0.5),
+        [1.0, 1.0, 1.0] => 0.0,
+        [-1.5, 1.0, -3.0] => sqrt(1.5^2 + 3^2),
+        [1.5, -1.0, 3.0] => sqrt(1.0^2 + 3^2),
+        [1.5, 1.0, -2.0] => 0.0,
+        [1.5, 1.0, -3.0] => 3 - sqrt(3) * sqrt(2),
+        [1.5, 1.0, 3.0] => 3 - sqrt(3) * sqrt(2);
+        mismatch = [1.0],
+    )
+    return
+end
+
+function test_normonecone()
+    _test_set(
+        MOI.NormOneCone(3),
+        [1.0, 1.0, 1.0] => 1.0,
+        [1.5, 1.0, -2.0] => 1.5,
+        [3.5, 1.0, -2.0] => 0.0;
+        mismatch = [1.0],
+    )
+    return
+end
+
+function test_norminfinitycone()
+    _test_set(
+        MOI.NormInfinityCone(3),
+        [1.0, 1.0, 1.0] => 0.0,
+        [1.5, 1.0, -2.0] => 0.5,
+        [3.5, 1.0, -2.0] => 0.0;
+        mismatch = [1.0],
+    )
+    return
+end
+
+function test_relativeentropycone()
+    _test_set(
+        MOI.RelativeEntropyCone(5),
+        [1.0, 1.0, 1.0, 1.0, 1.0] => 0.0,
+        [-2.0, 1.0, 1.0, 1.0, 1.0] => 2.0,
+        [1.0, 1.0, 2.0, 3.0, 1.0] => 3 * log(3 / 1) + 1 * log(1 / 2) - 1,
+        [4.0, 1.0, 2.0, 3.0, 1.0] => 0.0,
+        [4.0, -1.0, 2.0, 3.0, 1.0] => 2.0,
+        [4.0, 1.0, -2.0, 3.0, 1.0] => 3.0,
+        [4.0, 1.0, -2.0, 3.0, -1.0] => sqrt(3^2 + 2^2),
+        [0.0, 1.0, -2.0, 3.0, -1.0] => sqrt(3^2 + 2^2 + (3 * log(3))^2);
+        mismatch = [1.0],
+    )
+    return
+end
+
+function test_hyperrectangle()
+    _test_set(
+        MOI.HyperRectangle([0.0, 1.0], [1.0, 2.0]),
+        [0.0, 1.0] => 0.0,
+        [0.5, 1.2] => 0.0,
+        [-1.0, 1.5] => 1.0,
+        [0.5, 2.5] => 0.5,
+        [2.0, 0.0] => sqrt(2);
+        mismatch = [1.0],
+    )
+    return
+end
+
+function test_normcone()
+    _test_set(
+        MOI.NormCone(3, 4),
+        [1.0, 2.0, 3.0, 4.0] => LinearAlgebra.norm([2, 3, 4], 3) - 1,
+        [1.5, -2.0, 3.0, 4.0] => LinearAlgebra.norm([-2, 3, 4], 3) - 1.5;
+        mismatch = [1.0],
+    )
+    return
+end
+
+function test_sos1()
+    _test_set(
+        MOI.SOS1([1.0, 3.0, 2.0]),
+        [0.0, 1.0, 0.0] => 0.0,
+        [-0.5, 0.0, 0.0] => 0.0,
+        [1.0, 1.0, 0.0] => 1.0,
+        [-0.5, 1.5, 1.0] => sqrt(1 + 0.5^2);
+        mismatch = [1.0],
+    )
+    return
+end
+
+function test_sos2()
+    _test_set(
+        MOI.SOS2([1.0, 3.0, 2.0]),
+        [0.0, 1.0, 0.0] => 0.0,
+        [-0.5, 0.0, 0.0] => 0.0,
+        [0.0, 1.0, 1.0] => 0.0,
+        [-0.5, 0.0, 0.5] => 0.0,
+        [1.0, 1.0, 0.0] => 1.0,
+        [-0.5, 0.6, 0.0] => 0.5,
+        [-0.5, 1.5, 1.0] => 0.5;
+        mismatch = [1.0],
+    )
+    return
+end
+
 end
 
 TestFeasibilityChecker.runtests()