From 30c0a58cc3496e798d2b661f7a170d0af34caed4 Mon Sep 17 00:00:00 2001 From: odow Date: Tue, 24 Oct 2023 10:08:24 +1300 Subject: [PATCH] Add docstrings and better tests --- src/Utilities/distance_to_set.jl | 147 +++++++++++++++++++++++++----- test/Utilities/distance_to_set.jl | 113 +++++++++++------------ 2 files changed, 180 insertions(+), 80 deletions(-) diff --git a/src/Utilities/distance_to_set.jl b/src/Utilities/distance_to_set.jl index eabae3a569..d1848e9091 100644 --- a/src/Utilities/distance_to_set.jl +++ b/src/Utilities/distance_to_set.jl @@ -82,6 +82,7 @@ end ### MOI.AbstractScalarSets ### +# This is the minimal L2-norm. function distance_to_set( ::ProjectionUpperBoundDistance, x::T, @@ -90,6 +91,7 @@ function distance_to_set( return max(x - set.upper, zero(T)) end +# This is the minimal L2-norm. function distance_to_set( ::ProjectionUpperBoundDistance, x::T, @@ -98,6 +100,7 @@ function distance_to_set( return max(set.lower - x, zero(T)) end +# This is the minimal L2-norm. function distance_to_set( ::ProjectionUpperBoundDistance, x::T, @@ -106,6 +109,7 @@ function distance_to_set( return abs(set.value - x) end +# This is the minimal L2-norm. function distance_to_set( ::ProjectionUpperBoundDistance, x::T, @@ -114,6 +118,7 @@ function distance_to_set( return max(x - set.upper, set.lower - x, zero(T)) end +# This is the minimal L2-norm. function distance_to_set( ::ProjectionUpperBoundDistance, x::T, @@ -122,6 +127,7 @@ function distance_to_set( return min(abs(x - zero(T)), abs(x - one(T))) end +# This is the minimal L2-norm. function distance_to_set( ::ProjectionUpperBoundDistance, x::T, @@ -130,6 +136,7 @@ function distance_to_set( return abs(x - round(x)) end +# This is the minimal L2-norm. function distance_to_set( ::ProjectionUpperBoundDistance, x::T, @@ -138,6 +145,7 @@ function distance_to_set( return min(max(x - set.upper, set.lower - x, zero(T)), abs(x)) end +# This is the minimal L2-norm. function distance_to_set( ::ProjectionUpperBoundDistance, x::T, @@ -158,6 +166,7 @@ function _check_dimension(v::AbstractVector, s) return end +# This is the minimal L2-norm. function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, @@ -167,6 +176,7 @@ function distance_to_set( return LinearAlgebra.norm(max(-xi, zero(T)) for xi in x) end +# This is the minimal L2-norm. function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, @@ -176,6 +186,7 @@ function distance_to_set( return LinearAlgebra.norm(max(xi, zero(T)) for xi in x) end +# This is the minimal L2-norm. function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, @@ -185,6 +196,7 @@ function distance_to_set( return LinearAlgebra.norm(x) end +# This is the minimal L2-norm. function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, @@ -194,11 +206,13 @@ function distance_to_set( return zero(T) end +# This is the minimal L2-norm. function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, set::MOI.SecondOrderCone, ) where {T<:Real} + Base.require_one_based_indexing(x) _check_dimension(x, set) # Projections come from: # Parikh, N., & Boyd, S. (2014). Proximal algorithms. Foundations and @@ -214,26 +228,43 @@ function distance_to_set( return sqrt(2) / 2 * abs(t - rhs) end +""" + distance_to_set(::ProjectionUpperBoundDistance, x, ::MOI.RotatedSecondOrderCone) + +Let `(t, u, y...) = x`. Return the 2-norm of the vector `d` such that in `x + d`, +`u` is projected to `1` if `u <= 0`, and `t` is increased such that `x + d` +belongs to the set. +""" function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, set::MOI.RotatedSecondOrderCone, ) where {T<:Real} + Base.require_one_based_indexing(x) _check_dimension(x, set) - t, u, xs = x[1], x[2], @view(x[3:end]) + u, y = (x[2] <= 0 ? 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)), + u - x[2], # u -> 1 + max(LinearAlgebra.dot(y, y) / (2 * u) - x[1], zero(T)), ) return LinearAlgebra.norm(element_distance, 2) end +""" + distance_to_set(::ProjectionUpperBoundDistance, x, ::MOI.ExponentialCone) + +Let `(u, v, w) = x`. If `v > 0`, return the epigraph distance `d` such that +`(u, v, w + d)` belongs to the set. + +If `v <= 0` return the 2-norm of the vector `d` such that `x + d = (u, 1, z)` +where `z` satisfies the constraints. +""" function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, set::MOI.ExponentialCone, ) where {T<:Real} + Base.require_one_based_indexing(x) _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))) @@ -242,25 +273,44 @@ function distance_to_set( return max(x[2] * exp(x[1] / x[2]) - x[3], zero(T)) end +""" + distance_to_set(::ProjectionUpperBoundDistance, x, ::MOI.DualExponentialCone) + +Let `(u, v, w) = x`. If `u < 0`, return the epigraph distance `d` such that +`(u, v, w + d)` belongs to the set. + +If `u >= 0` return the 2-norm of the vector `d` such that `x + d = (u, -1, z)` +where `z` satisfies the constraints. +""" function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, set::MOI.DualExponentialCone, ) where {T<:Real} + Base.require_one_based_indexing(x) _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))) + element_distance = (x[1] - -one(T), max(exp(-x[2] - 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)) + return max(-x[1] * exp(x[2] / x[1] - 1) - x[3], zero(T)) end +""" + distance_to_set(::ProjectionUpperBoundDistance, x, ::MOI.GeometricMeanCone) + +Let `(t, y...) = x`. If all `y` are non-negative, return the epigraph distance +`d` such that `(t + d, y...)` belongs to the set. + +If any `y` are strictly negative, return the 2-norm of the vector `d` that +projects negative `y` elements to `0` and `t` to `ℝ₋`. +""" function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, set::MOI.GeometricMeanCone, ) where {T<:Real} + Base.require_one_based_indexing(x) _check_dimension(x, set) t, xs = x[1], @view(x[2:end]) if any(<(zero(T)), xs) # Project to x = 0 @@ -269,61 +319,105 @@ function distance_to_set( return max(t - prod(xs)^inv(MOI.dimension(set) - 1), zero(T)) end +""" + distance_to_set(::ProjectionUpperBoundDistance, x, ::MOI.PowerCone) + +Let `(a, b, c) = x`. If `a` and `b` are non-negative, return the epigraph +distance required to increase `c` such that the constraint is satisfied. + +If `a` or `b` is strictly negative, return the 2-norm of the vector `d` such +that in the vector `x + d`: `c`, and any negative `a` and `b` are projected to +`0`. +""" function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, set::MOI.PowerCone, ) where {T<:Real} + Base.require_one_based_indexing(x) _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, - ) + if x[1] >= 0 && x[2] >= 0 + return max(abs(x[3]) - x[1]^α * x[2]^(1 - α), zero(T)) end - return max(abs(x[3]) - x[1]^α * x[2]^(1 - α), zero(T)) + # Project to x = 0 + return LinearAlgebra.norm((min(x[1], zero(T)), min(x[2], zero(T)), x[3]), 2) end +""" + distance_to_set(::ProjectionUpperBoundDistance, x, ::MOI.DualPowerCone) + +Let `(a, b, c) = x`. If `a` and `b` are non-negative, return the epigraph +distance required to increase `c` such that the constraint is satisfied. + +If `a` or `b` is strictly negative, return the 2-norm of the vector `d` such +that in the vector `x + d`: `c`, and any negative `a` and `b` are projected to +`0`. +""" function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, set::MOI.DualPowerCone, ) where {T<:Real} + Base.require_one_based_indexing(x) _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, - ) + if x[1] >= 0 && x[2] >= 0 + return max(abs(x[3]) - (x[1] / α)^α * (x[2] / (1 - α))^(1 - α), zero(T)) end - return max(abs(x[3]) - (x[1] / α)^α * (x[2] / (1 - α))^(1 - α), zero(T)) + # Project to x = 0 + return LinearAlgebra.norm((min(x[1], zero(T)), min(x[2], zero(T)), x[3]), 2) end +""" + distance_to_set(::ProjectionUpperBoundDistance, x, ::MOI.NormOneCone) + +Let `(t, y...) = x`. Return the epigraph distance `d` such that `(t + d, y...)` +belongs to the set. +""" function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, set::MOI.NormOneCone, ) where {T<:Real} + Base.require_one_based_indexing(x) _check_dimension(x, set) return max(sum(abs, @view(x[2:end])) - x[1], zero(T)) end +""" + distance_to_set(::ProjectionUpperBoundDistance, x, ::MOI.NormInfinityCone) + +Let `(t, y...) = x`. Return the epigraph distance `d` such that `(t + d, y...)` +belongs to the set. +""" function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, set::MOI.NormInfinityCone, ) where {T<:Real} + Base.require_one_based_indexing(x) _check_dimension(x, set) return max(maximum(abs, @view(x[2:end])) - x[1], zero(T)) end +""" + distance_to_set(::ProjectionUpperBoundDistance, x, ::MOI.RelativeEntropyCone) + +Let `(u, v..., w...) = x`. If `v` and `w` are strictly positive, return the +epigraph distance required to increase `u` such that the constraint is satisfied. + +If any elements in `v` or `w` are non-positive, return the 2-norm of the vector +`d` such that in the vector `x + d`: any non-positive elements in `v` and `w` +are projected to `1`, and `u` is projected such that the epigraph constraint +holds. +""" function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, set::MOI.RelativeEntropyCone, ) where {T<:Real} + Base.require_one_based_indexing(x) _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]) @@ -343,6 +437,7 @@ function distance_to_set( return max(sum(w[i] * log(w[i] / v[i]) for i in eachindex(w)) - u, zero(T)) end +# This is the minimal L2-norm. function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, @@ -354,15 +449,23 @@ function distance_to_set( return LinearAlgebra.norm(element_distance, 2) end +""" + distance_to_set(::ProjectionUpperBoundDistance, x, set::MOI.NormCone) + +Let `(t, y...) = x`. Return the epigraph distance `d` such that `(t + d, y...)` +belongs to the set. +""" function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, set::MOI.NormCone, ) where {T<:Real} + Base.require_one_based_indexing(x) _check_dimension(x, set) return max(LinearAlgebra.norm(@view(x[2:end]), set.p) - x[1], zero(T)) end +# This is the minimal L2-norm. function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, @@ -370,9 +473,10 @@ function distance_to_set( ) where {T<:Real} _check_dimension(x, set) _, i = findmax(abs, x) - return LinearAlgebra.norm2([x[j] for j in eachindex(x) if j != i]) + return LinearAlgebra.norm([x[j] for j in eachindex(x) if j != i], 2) end +# This is the minimal L2-norm. function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, @@ -382,5 +486,6 @@ function distance_to_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])]) + elements = [x[i] for i in eachindex(x) if !(i in pairs[k])] + return LinearAlgebra.norm(elements, 2) end diff --git a/test/Utilities/distance_to_set.jl b/test/Utilities/distance_to_set.jl index 7ef315f146..558080efaf 100644 --- a/test/Utilities/distance_to_set.jl +++ b/test/Utilities/distance_to_set.jl @@ -44,113 +44,107 @@ function test_unsupported() end function test_lessthan() - @test MOI.Utilities.distance_to_set(1.0, MOI.LessThan(2.0)) ≈ 0.0 - @test MOI.Utilities.distance_to_set(1.0, MOI.LessThan(0.5)) ≈ 0.5 + _test_set(MOI.LessThan(2.0), 1.0 => 0.0) + _test_set(MOI.LessThan(0.5), 1.0 => 0.5) return end function test_greaterthan() - @test MOI.Utilities.distance_to_set(1.0, MOI.GreaterThan(2.0)) ≈ 1.0 - @test MOI.Utilities.distance_to_set(1.0, MOI.GreaterThan(0.5)) ≈ 0.0 + _test_set(MOI.GreaterThan(2.0), 1.0 => 1.0) + _test_set(MOI.GreaterThan(0.5), 1.0 => 0.0) return end function test_equalto() - @test MOI.Utilities.distance_to_set(1.0, MOI.EqualTo(2.0)) ≈ 1.0 - @test MOI.Utilities.distance_to_set(1.0, MOI.EqualTo(0.5)) ≈ 0.5 + _test_set(MOI.EqualTo(2.0), 1.0 => 1.0) + _test_set(MOI.EqualTo(0.5), 1.0 => 0.5) return end function test_interval() - @test MOI.Utilities.distance_to_set(1.0, MOI.Interval(1.0, 2.0)) ≈ 0.0 - @test MOI.Utilities.distance_to_set(0.5, MOI.Interval(1.0, 2.0)) ≈ 0.5 - @test MOI.Utilities.distance_to_set(2.75, MOI.Interval(1.0, 2.0)) ≈ 0.75 + _test_set(MOI.Interval(1.0, 2.0), 1.0 => 0.0, 0.5 => 0.5, 2.75 => 0.75) return end function test_zeroone() - @test MOI.Utilities.distance_to_set(0.6, MOI.ZeroOne()) ≈ 0.4 - @test MOI.Utilities.distance_to_set(-0.01, MOI.ZeroOne()) ≈ 0.01 - @test MOI.Utilities.distance_to_set(1.01, MOI.ZeroOne()) ≈ 0.01 + _test_set( + MOI.ZeroOne(), + 0.0 => 0.0, + 1.0 => 0.0, + 2.0 => 1.0, + 0.6 => 0.4, + -0.01 => 0.01, + 1.01 => 0.01, + ) return end function test_integer() - @test MOI.Utilities.distance_to_set(0.6, MOI.Integer()) ≈ 0.4 - @test MOI.Utilities.distance_to_set(3.1, MOI.Integer()) ≈ 0.1 - @test MOI.Utilities.distance_to_set(-0.01, MOI.Integer()) ≈ 0.01 - @test MOI.Utilities.distance_to_set(1.01, MOI.Integer()) ≈ 0.01 + _test_set( + MOI.Integer(), + 0.6 => 0.4, + 3.1 => 0.1, + -0.01 => 0.01, + 1.01 => 0.01, + ) return end function test_semicontinuous() - s = MOI.Semicontinuous(2.0, 4.0) - @test MOI.Utilities.distance_to_set(-2.0, s) ≈ 2.0 - @test MOI.Utilities.distance_to_set(0.5, s) ≈ 0.5 - @test MOI.Utilities.distance_to_set(1.9, s) ≈ 0.1 - @test MOI.Utilities.distance_to_set(2.1, s) ≈ 0.0 - @test MOI.Utilities.distance_to_set(4.1, s) ≈ 0.1 + _test_set( + MOI.Semicontinuous(2.0, 4.0), + -2.0 => 2.0, + 0.5 => 0.5, + 1.9 => 0.1, + 2.1 => 0.0, + 4.1 => 0.1, + ) return end function test_semiintger() - s = MOI.Semiinteger(1.9, 4.0) - @test MOI.Utilities.distance_to_set(-2.0, s) ≈ 2.0 - @test MOI.Utilities.distance_to_set(0.5, s) ≈ 0.5 - @test MOI.Utilities.distance_to_set(1.9, s) ≈ 0.1 - @test MOI.Utilities.distance_to_set(2.1, s) ≈ 0.1 - @test MOI.Utilities.distance_to_set(4.1, s) ≈ 0.1 + _test_set( + MOI.Semiinteger(1.9, 4.0), + -2.0 => 2.0, + 0.5 => 0.5, + 1.9 => 0.1, + 2.1 => 0.1, + 4.1 => 0.1, + ) return end function test_nonnegatives() - @test_throws( - DimensionMismatch, - MOI.Utilities.distance_to_set([-1.0, 1.0], MOI.Nonnegatives(1)) - ) - @test MOI.Utilities.distance_to_set([-1.0, 1.0], MOI.Nonnegatives(2)) ≈ 1.0 + _test_set(MOI.Nonnegatives(2), [-1.0, 1.0] => 1.0; mismatch = [1.0]) return end function test_nonpositives() - @test_throws( - DimensionMismatch, - MOI.Utilities.distance_to_set([-1.0, 1.0], MOI.Nonpositives(1)) - ) - @test MOI.Utilities.distance_to_set([-1.0, 1.0], MOI.Nonpositives(2)) ≈ 1.0 + _test_set(MOI.Nonpositives(2), [-1.0, 1.0] => 1.0; mismatch = [1.0]) return end function test_reals() - @test_throws( - DimensionMismatch, - MOI.Utilities.distance_to_set([-1.0, 1.0], MOI.Reals(1)) - ) - @test MOI.Utilities.distance_to_set([-1.0, 1.0], MOI.Reals(2)) ≈ 0.0 + _test_set(MOI.Reals(2), [-1.0, 1.0] => 0.0; mismatch = [1.0]) return end function test_zeros() - @test_throws( - DimensionMismatch, - MOI.Utilities.distance_to_set([-1.0, 1.0], MOI.Zeros(1)) - ) - @test MOI.Utilities.distance_to_set([-1.0, 1.0], MOI.Zeros(2)) ≈ sqrt(2) + _test_set(MOI.Zeros(2), [-1.0, 1.0] => sqrt(2); mismatch = [1.0]) return end function test_secondordercone() - @test_throws( - DimensionMismatch, - MOI.Utilities.distance_to_set([-1.0, 1.0], MOI.SecondOrderCone(3)) - ) - set = MOI.SecondOrderCone(3) - @test MOI.Utilities.distance_to_set([sqrt(2), 1, 1], set) ≈ 0.0 - @test MOI.Utilities.distance_to_set([-sqrt(2), 1, 1], set) ≈ 2 - @test MOI.Utilities.distance_to_set([-2, 1, 1], set) ≈ sqrt(6) # According to Boyd, (t, x) = (1, [1, 1]), projects to: d = ((1 / 2) * (1 + 1 / √2) * [√2, 1, 1]) .- [1, 1, 1] - @test MOI.Utilities.distance_to_set([1, 1, 1], set) ≈ LinearAlgebra.norm(d) + _test_set( + MOI.SecondOrderCone(3), + [sqrt(2), 1, 1] => 0.0, + [-sqrt(2), 1, 1] => 2, + [-2, 1, 1] => sqrt(6), + [1, 1, 1] => LinearAlgebra.norm(d); + mismatch = [-1.0, 1.0], + ) return end @@ -158,7 +152,8 @@ 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); + [-1.0, 1.0, 1.0, 1.0] => 2.0, + [1.0, 0.0, 2.0, 3.0] => sqrt(1 + 5.5^2); mismatch = [1.0], ) return @@ -180,7 +175,7 @@ function test_dualexponential() 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); + [-2.0, 3.0, 0.1] => 2 * exp(3 / -2 - 1) - 0.1; mismatch = [1.0], ) return