From b606dc3dda40db5f49b5309d79a43fc2a48938ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Tue, 11 Feb 2020 10:55:17 +0100 Subject: [PATCH 01/28] first sets --- src/Utilities/Utilities.jl | 1 + src/Utilities/set_eval.jl | 52 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 53 insertions(+) create mode 100644 src/Utilities/set_eval.jl diff --git a/src/Utilities/Utilities.jl b/src/Utilities/Utilities.jl index 53c5d3477b..1a64e077e5 100644 --- a/src/Utilities/Utilities.jl +++ b/src/Utilities/Utilities.jl @@ -32,6 +32,7 @@ include("constraints.jl") include("copy.jl") include("results.jl") include("variables.jl") +include("set_eval.jl") include("model.jl") include("parser.jl") diff --git a/src/Utilities/set_eval.jl b/src/Utilities/set_eval.jl new file mode 100644 index 0000000000..73c734411c --- /dev/null +++ b/src/Utilities/set_eval.jl @@ -0,0 +1,52 @@ + +Base.in(v, s::MOI.LessThan) = v <= s.upper +Base.in(v, s::MOI.GreaterThan) = v >= s.lower +Base.in(v, s::MOI.EqualTo) = v ≈ s.value + +Base.in(v, s::MOI.Interval) = s.lower <= v <= s.upper + +Base.in(v::AbstractVector{<:Real}, s::MOI.Reals) = _dim_match(v, s) +Base.in(v, s::MOI.Zeros) = _dim_match(v, s) && all(≈(0), v) + +Base.in(v::AbstractVector{<:Real}, s::MOI.Nonnegatives) = _dim_match(v, s) && all(≥(0), v) +Base.in(v::AbstractVector{<:Real}, s::MOI.Nonpositives) = _dim_match(v, s) && all(≤(0), v) + +# Norm cones + +function Base.in(v::AbstractVector{<:Real}, s::MOI.NormInfinityCone) + _dim_match(v, s) || return false + t = first(v) + xs = v[2:end] + for x in xs + if abs(x) > t + return false + end + end + return true +end + +function Base.in(v::AbstractVector{<:Real}, s::MOI.NormOneCone) + _dim_match(v, s) || return false + t = v[1] + xs = v[2:end] + return t >= sum(abs, xs) +end + +function Base.in(v::AbstractVector{<:Real}, s::MOI.SecondOrderCone) + _dim_match(v, s) || return false + t = v[1] + xs = v[2:end] + return t ≥ 0 && t^2 >= dot(xs, xs) # avoids sqrt +end + + +function Base.in(v::AbstractVector{<:Real}, s::MOI.RotatedSecondOrderCone) + _dim_match(v, s) || return false + t = v[1] + u = v[2] + t ≥ 0 && u ≥ 0 || return false + xs = v[3:end] + return 2 * t * u >= dot(xs, xs) +end + +_dim_match(v::AbstractVector, s::MOI.AbstractVectorSet) = length(v) == MOI.dimension(s) From 6670920f3f7f25d881d4a56b05b0b63b2f60b4bb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Tue, 11 Feb 2020 11:40:28 +0100 Subject: [PATCH 02/28] more sets --- src/Utilities/set_eval.jl | 93 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 92 insertions(+), 1 deletion(-) diff --git a/src/Utilities/set_eval.jl b/src/Utilities/set_eval.jl index 73c734411c..70c785d8fc 100644 --- a/src/Utilities/set_eval.jl +++ b/src/Utilities/set_eval.jl @@ -5,6 +5,8 @@ Base.in(v, s::MOI.EqualTo) = v ≈ s.value Base.in(v, s::MOI.Interval) = s.lower <= v <= s.upper +_dim_match(v::AbstractVector, s::MOI.AbstractVectorSet) = length(v) == MOI.dimension(s) + Base.in(v::AbstractVector{<:Real}, s::MOI.Reals) = _dim_match(v, s) Base.in(v, s::MOI.Zeros) = _dim_match(v, s) && all(≈(0), v) @@ -49,4 +51,93 @@ function Base.in(v::AbstractVector{<:Real}, s::MOI.RotatedSecondOrderCone) return 2 * t * u >= dot(xs, xs) end -_dim_match(v::AbstractVector, s::MOI.AbstractVectorSet) = length(v) == MOI.dimension(s) +function Base.in(v::AbstractVector{<:Real}, s::MOI.GeometricMeanCone) + t = v[1] + xs = v[2:end] + n = MOI.dimension(s) - 1 + xs in MOI.Nonnegatives(n) || return false # this also checks dimensions + return t^n <= prod(xs) +end + +function Base.in(v::AbstractVector{<:Real}, s::MOI.ExponentialCone) + length(v) == 3 || return false + @inbounds y = v[2] + y > 0 || return false + @inbounds x = v[1] + @inbounds z = v[3] + return y * exp(x/y) <= z +end + +function Base.in(v::AbstractVector{<:Real}, s::MOI.DualExponentialCone) + length(v) == 3 || return false + @inbounds u = v[1] + u < 0 || return false + @inbounds v = v[2] + @inbounds w = v[3] + return -u*exp(v/u) ≤ ℯ * w +end + +function Base.in(v::AbstractVector{<:Real}, s::MOI.PowerCone) + length(v) == 3 || return false + @inbounds x = v[1] + x ≥ 0 || return false + @inbounds y = v[2] + y ≥ 0 || return false + @inbounds z = v[3] + e = s.exponent + return x^e * y^(1-e) ≥ abs(z) +end + +function Base.in(v::AbstractVector{<:Real}, s::MOI.DualPowerCone) + length(v) == 3 || return false + @inbounds u = v[1] + u ≥ 0 || return false + @inbounds v = v[2] + v ≥ 0 || return false + @inbounds w = v[3] + e = s.exponent + ce = 1-e + return (u/e)^e * (v/ce)^ce ≥ abs(w) +end + +function Base.in(v::AbstractVector{<:Real}, s::MOI.RelativeEntropyCone) + _dim_match(v, s) || return false + all(>=(0), v[2:end]) || return false + n = (MOI.dimension(s)-1) ÷ 2 + @inbounds u = v[1] + v = [2:n+1] + w = [n+2:end] + @inbounds s = sum(w[i] * log(w[i]/v[i]) for i in eachindex(w)) + return u ≥ s +end + + +function Base.in(v::AbstractVector{<:Real}, s::MOI.NormSpectralCone) + _dim_match(v, s) || return false + t = v[1] + m = reshape(v[2:end], (s.row_dim, s.column_dim)) + s1 = LinearAlgebra.svd(m).S[1] + return t ≥ s1 +end + +function Base.in(v::AbstractVector{<:Real}, s::MOI.NormNuclearCone) + _dim_match(v, s) || return false + t = v[1] + m = reshape(v[2:end], (s.row_dim, s.column_dim)) + s = sum(LinearAlgebra.svd(m).S) + return t ≥ s +end + +function Base.in(v::AbstractVector{<:Real}, s::MOI.SOS1) + _dim_match(v, s) || return false + found_nonzero = false + for x in v + if !isapprox(x, 0) + if found_nonzero + return false + end + found_nonzero = true + end + end + return true +end From 93ac7f9d96b1906fc9e3130fe8747dac791fd8ed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Tue, 11 Feb 2020 11:48:31 +0100 Subject: [PATCH 03/28] forgot subslice --- src/Utilities/set_eval.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Utilities/set_eval.jl b/src/Utilities/set_eval.jl index 70c785d8fc..6fb8a621b6 100644 --- a/src/Utilities/set_eval.jl +++ b/src/Utilities/set_eval.jl @@ -105,8 +105,8 @@ function Base.in(v::AbstractVector{<:Real}, s::MOI.RelativeEntropyCone) all(>=(0), v[2:end]) || return false n = (MOI.dimension(s)-1) ÷ 2 @inbounds u = v[1] - v = [2:n+1] - w = [n+2:end] + v = v[2:(n+1)] + w = v[(n+2):end] @inbounds s = sum(w[i] * log(w[i]/v[i]) for i in eachindex(w)) return u ≥ s end From 55ed9d6b5cf45179ac4ed4c96e706ec89c1cd8f9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Tue, 11 Feb 2020 17:04:39 +0100 Subject: [PATCH 04/28] rewrite as distance function --- src/Utilities/set_eval.jl | 139 ++++++++++++++++++++------------------ 1 file changed, 73 insertions(+), 66 deletions(-) diff --git a/src/Utilities/set_eval.jl b/src/Utilities/set_eval.jl index 6fb8a621b6..6da0e5253d 100644 --- a/src/Utilities/set_eval.jl +++ b/src/Utilities/set_eval.jl @@ -1,114 +1,121 @@ -Base.in(v, s::MOI.LessThan) = v <= s.upper -Base.in(v, s::MOI.GreaterThan) = v >= s.lower -Base.in(v, s::MOI.EqualTo) = v ≈ s.value +function set_distance(v::Real, s) + d = _distance(v, s) + return d < 0 : zero(d) : d +end + +function set_distance(v::AbstractVector{T}, s) where {T <: Real} + length(v) != MOI.dimension(s) && throw(DimensionMismatch("Mismatch between v and s")) + d = _distance(v, s) + return [di < 0 : zero(T) : di for di in d] +end + +_distance(v, s::MOI.LessThan) = s.upper - v +_distance(v, s::MOI.GreaterThan) = v - s.lower +_distance(v, s::MOI.EqualTo) = abs(v - s.value) -Base.in(v, s::MOI.Interval) = s.lower <= v <= s.upper +_distance(v, s::MOI.Interval) = max(v - s.lower, s.upper - v) -_dim_match(v::AbstractVector, s::MOI.AbstractVectorSet) = length(v) == MOI.dimension(s) +_distance(v::AbstractVector{<:Real}, ::MOI.Reals) = false * v -Base.in(v::AbstractVector{<:Real}, s::MOI.Reals) = _dim_match(v, s) -Base.in(v, s::MOI.Zeros) = _dim_match(v, s) && all(≈(0), v) +_distance(v, ::MOI.Zeros) = abs.(v) -Base.in(v::AbstractVector{<:Real}, s::MOI.Nonnegatives) = _dim_match(v, s) && all(≥(0), v) -Base.in(v::AbstractVector{<:Real}, s::MOI.Nonpositives) = _dim_match(v, s) && all(≤(0), v) +_distance(v::AbstractVector{<:Real}, ::MOI.Nonnegatives) = -v +_distance(v::AbstractVector{<:Real}, ::MOI.Nonpositives) = v # Norm cones -function Base.in(v::AbstractVector{<:Real}, s::MOI.NormInfinityCone) - _dim_match(v, s) || return false +function _distance(v::AbstractVector{<:Real}, ::MOI.NormInfinityCone) t = first(v) xs = v[2:end] - for x in xs - if abs(x) > t - return false - end - end - return true + return [t - abs(x) for x in xs] end -function Base.in(v::AbstractVector{<:Real}, s::MOI.NormOneCone) - _dim_match(v, s) || return false +function _distance(v::AbstractVector{<:Real}, ::MOI.NormOneCone) t = v[1] xs = v[2:end] - return t >= sum(abs, xs) + return t - sum(abs, xs) end -function Base.in(v::AbstractVector{<:Real}, s::MOI.SecondOrderCone) - _dim_match(v, s) || return false +function _distance(v::AbstractVector{<:Real}, ::MOI.SecondOrderCone) t = v[1] xs = v[2:end] - return t ≥ 0 && t^2 >= dot(xs, xs) # avoids sqrt + return min(-t, dot(xs, xs) - t^2) # avoids sqrt end - -function Base.in(v::AbstractVector{<:Real}, s::MOI.RotatedSecondOrderCone) - _dim_match(v, s) || return false +function _distance(v::AbstractVector{<:Real}, ::MOI.RotatedSecondOrderCone) t = v[1] u = v[2] - t ≥ 0 && u ≥ 0 || return false xs = v[3:end] - return 2 * t * u >= dot(xs, xs) + return min( + -t, -u, + dot(xs, xs) - 2 * t * u + ) end -function Base.in(v::AbstractVector{<:Real}, s::MOI.GeometricMeanCone) +function _distance(v::AbstractVector{<:Real}, s::MOI.GeometricMeanCone) t = v[1] xs = v[2:end] n = MOI.dimension(s) - 1 - xs in MOI.Nonnegatives(n) || return false # this also checks dimensions - return t^n <= prod(xs) + return min( + minimum(xs), + t^n - prod(xs), + ) end -function Base.in(v::AbstractVector{<:Real}, s::MOI.ExponentialCone) - length(v) == 3 || return false - @inbounds y = v[2] - y > 0 || return false - @inbounds x = v[1] - @inbounds z = v[3] - return y * exp(x/y) <= z +function _distance(v::AbstractVector{<:Real}, ::MOI.ExponentialCone) + x = v[1] + y = v[2] + z = v[3] + return min( + y, + y * exp(x/y) - z, + ) end -function Base.in(v::AbstractVector{<:Real}, s::MOI.DualExponentialCone) - length(v) == 3 || return false - @inbounds u = v[1] - u < 0 || return false - @inbounds v = v[2] - @inbounds w = v[3] - return -u*exp(v/u) ≤ ℯ * w +function _distance(v::AbstractVector{<:Real}, ::MOI.DualExponentialCone) + u = v[1] + v = v[2] + w = v[3] + return min( + -u, + -u*exp(v/u) - ℯ * w + ) end -function Base.in(v::AbstractVector{<:Real}, s::MOI.PowerCone) - length(v) == 3 || return false - @inbounds x = v[1] - x ≥ 0 || return false - @inbounds y = v[2] - y ≥ 0 || return false - @inbounds z = v[3] +function _distance(v::AbstractVector{<:Real}, s::MOI.PowerCone) + x = v[1] + y = v[2] + z = v[3] e = s.exponent - return x^e * y^(1-e) ≥ abs(z) + return min( + x, + y, + abs(z) - x^e * y^(1-e) + ) end -function Base.in(v::AbstractVector{<:Real}, s::MOI.DualPowerCone) - length(v) == 3 || return false - @inbounds u = v[1] - u ≥ 0 || return false - @inbounds v = v[2] - v ≥ 0 || return false - @inbounds w = v[3] +function _distance(v::AbstractVector{<:Real}, s::MOI.DualPowerCone) + u = v[1] + v = v[2] + w = v[3] e = s.exponent ce = 1-e - return (u/e)^e * (v/ce)^ce ≥ abs(w) + return min( + u, + v, + abs(w) - (u/e)^e * (v/ce)^ce + ) end -function Base.in(v::AbstractVector{<:Real}, s::MOI.RelativeEntropyCone) - _dim_match(v, s) || return false +function _distance(v::AbstractVector{<:Real}, s::MOI.RelativeEntropyCone) all(>=(0), v[2:end]) || return false n = (MOI.dimension(s)-1) ÷ 2 - @inbounds u = v[1] + u = v[1] v = v[2:(n+1)] w = v[(n+2):end] - @inbounds s = sum(w[i] * log(w[i]/v[i]) for i in eachindex(w)) - return u ≥ s + s = sum(w[i] * log(w[i]/v[i]) for i in eachindex(w)) + return s - u end From 09663ddd88352516fba9515b3c3bc0ce326082c8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Tue, 11 Feb 2020 17:13:29 +0100 Subject: [PATCH 05/28] sign confusion --- src/Utilities/set_eval.jl | 34 ++++++++++++++++++---------------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/src/Utilities/set_eval.jl b/src/Utilities/set_eval.jl index 6da0e5253d..4bcd9a62a1 100644 --- a/src/Utilities/set_eval.jl +++ b/src/Utilities/set_eval.jl @@ -10,11 +10,11 @@ function set_distance(v::AbstractVector{T}, s) where {T <: Real} return [di < 0 : zero(T) : di for di in d] end -_distance(v, s::MOI.LessThan) = s.upper - v -_distance(v, s::MOI.GreaterThan) = v - s.lower +_distance(v, s::MOI.LessThan) = v - s.upper +_distance(v, s::MOI.GreaterThan) = s.lower - v _distance(v, s::MOI.EqualTo) = abs(v - s.value) -_distance(v, s::MOI.Interval) = max(v - s.lower, s.upper - v) +_distance(v, s::MOI.Interval) = max(s.lower - v, v - s.upper) _distance(v::AbstractVector{<:Real}, ::MOI.Reals) = false * v @@ -40,14 +40,14 @@ end function _distance(v::AbstractVector{<:Real}, ::MOI.SecondOrderCone) t = v[1] xs = v[2:end] - return min(-t, dot(xs, xs) - t^2) # avoids sqrt + return max(-t, dot(xs, xs) - t^2) # avoids sqrt end function _distance(v::AbstractVector{<:Real}, ::MOI.RotatedSecondOrderCone) t = v[1] u = v[2] xs = v[3:end] - return min( + return max( -t, -u, dot(xs, xs) - 2 * t * u ) @@ -57,8 +57,8 @@ function _distance(v::AbstractVector{<:Real}, s::MOI.GeometricMeanCone) t = v[1] xs = v[2:end] n = MOI.dimension(s) - 1 - return min( - minimum(xs), + return max( + -minimum(xs), t^n - prod(xs), ) end @@ -67,7 +67,7 @@ function _distance(v::AbstractVector{<:Real}, ::MOI.ExponentialCone) x = v[1] y = v[2] z = v[3] - return min( + return max( y, y * exp(x/y) - z, ) @@ -77,7 +77,7 @@ function _distance(v::AbstractVector{<:Real}, ::MOI.DualExponentialCone) u = v[1] v = v[2] w = v[3] - return min( + return max( -u, -u*exp(v/u) - ℯ * w ) @@ -88,7 +88,7 @@ function _distance(v::AbstractVector{<:Real}, s::MOI.PowerCone) y = v[2] z = v[3] e = s.exponent - return min( + return max( x, y, abs(z) - x^e * y^(1-e) @@ -101,26 +101,28 @@ function _distance(v::AbstractVector{<:Real}, s::MOI.DualPowerCone) w = v[3] e = s.exponent ce = 1-e - return min( + return max( u, v, abs(w) - (u/e)^e * (v/ce)^ce ) end -function _distance(v::AbstractVector{<:Real}, s::MOI.RelativeEntropyCone) +function _distance(v::AbstractVector{<:Real}, set::MOI.RelativeEntropyCone) all(>=(0), v[2:end]) || return false - n = (MOI.dimension(s)-1) ÷ 2 + n = (MOI.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)) - return s - u + return max( + -minimum(v[2:end]), + s - u, + ) end -function Base.in(v::AbstractVector{<:Real}, s::MOI.NormSpectralCone) - _dim_match(v, s) || return false +function _distance(v::AbstractVector{<:Real}, s::MOI.NormSpectralCone) t = v[1] m = reshape(v[2:end], (s.row_dim, s.column_dim)) s1 = LinearAlgebra.svd(m).S[1] From 9d4def8cde3154e7e859c7165c3c2ba853784a7e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Tue, 11 Feb 2020 17:53:03 +0100 Subject: [PATCH 06/28] added integers --- src/Utilities/set_eval.jl | 53 ++++++++++++++++++++++++++++----------- 1 file changed, 38 insertions(+), 15 deletions(-) diff --git a/src/Utilities/set_eval.jl b/src/Utilities/set_eval.jl index 4bcd9a62a1..3742b60bbe 100644 --- a/src/Utilities/set_eval.jl +++ b/src/Utilities/set_eval.jl @@ -1,10 +1,10 @@ -function set_distance(v::Real, s) +function distance_to_set(v::Real, s) d = _distance(v, s) return d < 0 : zero(d) : d end -function set_distance(v::AbstractVector{T}, s) where {T <: Real} +function distance_to_set(v::AbstractVector{T}, s) where {T <: Real} length(v) != MOI.dimension(s) && throw(DimensionMismatch("Mismatch between v and s")) d = _distance(v, s) return [di < 0 : zero(T) : di for di in d] @@ -126,27 +126,50 @@ function _distance(v::AbstractVector{<:Real}, s::MOI.NormSpectralCone) t = v[1] m = reshape(v[2:end], (s.row_dim, s.column_dim)) s1 = LinearAlgebra.svd(m).S[1] - return t ≥ s1 + return s1 - t end -function Base.in(v::AbstractVector{<:Real}, s::MOI.NormNuclearCone) - _dim_match(v, s) || return false +function _distance(v::AbstractVector{<:Real}, s::MOI.NormNuclearCone) t = v[1] m = reshape(v[2:end], (s.row_dim, s.column_dim)) - s = sum(LinearAlgebra.svd(m).S) - return t ≥ s + s1 = sum(LinearAlgebra.svd(m).S) + return s1 - t end -function Base.in(v::AbstractVector{<:Real}, s::MOI.SOS1) - _dim_match(v, s) || return false - found_nonzero = false - for x in v +# return the second largest absolute value (=0 if SOS1 respected) +function _distance(v::AbstractVector{T}, s::MOI.SOS1) where {T <: Real} + first_non_zero = second_non_zero = -1 + v1 = zero(T) + v2 = zero(T) + for (i, x) in enumerate(v) if !isapprox(x, 0) - if found_nonzero - return false + if first_non_zero < 0 + first_non_zero = i + v1 = x + elseif second_non_zero < 0 + second_non_zero = i + v2 = x + else # replace one of the two + xa = abs(x) + if xa > v1 + second_non_zero = first_non_zero + v2 = v1 + first_non_zero = i + v1 = xa + elseif v[second_non_zero] <= x < v[first_non_zero] + second_non_zero = i + v2 = xa + end end - found_nonzero = true end end - return true + return v2 +end + +function _distance(v::T, ::MOI.ZeroOne) where {T <: Real} + return min(abs(v - zero(T)), abs(v - one(T))) +end + +function _distance(v::Real, ::MOI.Integer) + return min(abs(v - floor(v)), abs(v - ceil(v))) end From c02e1badf35446d3f14c1aa06f2a94a098abe902 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Tue, 11 Feb 2020 23:23:00 +0100 Subject: [PATCH 07/28] function for verification not being toplevel --- src/Utilities/set_eval.jl | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/src/Utilities/set_eval.jl b/src/Utilities/set_eval.jl index 3742b60bbe..0c054b16ca 100644 --- a/src/Utilities/set_eval.jl +++ b/src/Utilities/set_eval.jl @@ -1,24 +1,27 @@ -function distance_to_set(v::Real, s) - d = _distance(v, s) - return d < 0 : zero(d) : d +function _correct_distance_to_set(d::Real) + return ifelse(d < 0, zero(d), d) +end + +function _correct_distance_to_set(d::AbstractVector{T}) where {T <: Real} + return [ifelse(di < 0, zero(T), di) for di in d] end -function distance_to_set(v::AbstractVector{T}, s) where {T <: Real} +function distance_to_set(v, s) length(v) != MOI.dimension(s) && throw(DimensionMismatch("Mismatch between v and s")) d = _distance(v, s) - return [di < 0 : zero(T) : di for di in d] + return _correct_distance_to_set(d) end -_distance(v, s::MOI.LessThan) = v - s.upper -_distance(v, s::MOI.GreaterThan) = s.lower - v -_distance(v, s::MOI.EqualTo) = abs(v - s.value) +_distance(v::Real, s::MOI.LessThan) = v - s.upper +_distance(v::Real, s::MOI.GreaterThan) = s.lower - v +_distance(v::Real, s::MOI.EqualTo) = abs(v - s.value) -_distance(v, s::MOI.Interval) = max(s.lower - v, v - s.upper) +_distance(v::Real, s::MOI.Interval) = max(s.lower - v, v - s.upper) _distance(v::AbstractVector{<:Real}, ::MOI.Reals) = false * v -_distance(v, ::MOI.Zeros) = abs.(v) +_distance(v::AbstractVector{<:Real}, ::MOI.Zeros) = abs.(v) _distance(v::AbstractVector{<:Real}, ::MOI.Nonnegatives) = -v _distance(v::AbstractVector{<:Real}, ::MOI.Nonpositives) = v From eea190e551434cdd62c683916322b0a853704f87 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Wed, 12 Feb 2020 09:26:52 +0100 Subject: [PATCH 08/28] more integer sets --- src/Utilities/set_eval.jl | 26 +++++++++++++++++++++----- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/src/Utilities/set_eval.jl b/src/Utilities/set_eval.jl index 0c054b16ca..097655d66b 100644 --- a/src/Utilities/set_eval.jl +++ b/src/Utilities/set_eval.jl @@ -139,6 +139,16 @@ function _distance(v::AbstractVector{<:Real}, s::MOI.NormNuclearCone) return s1 - t end +## Integer sets + +function _distance(v::T, ::MOI.ZeroOne) where {T <: Real} + return min(abs(v - zero(T)), abs(v - one(T))) +end + +function _distance(v::Real, ::MOI.Integer) + return min(abs(v - floor(v)), abs(v - ceil(v))) +end + # return the second largest absolute value (=0 if SOS1 respected) function _distance(v::AbstractVector{T}, s::MOI.SOS1) where {T <: Real} first_non_zero = second_non_zero = -1 @@ -169,10 +179,16 @@ function _distance(v::AbstractVector{T}, s::MOI.SOS1) where {T <: Real} return v2 end -function _distance(v::T, ::MOI.ZeroOne) where {T <: Real} - return min(abs(v - zero(T)), abs(v - one(T))) -end +# TODO SOS2 -function _distance(v::Real, ::MOI.Integer) - return min(abs(v - floor(v)), abs(v - ceil(v))) +# takes in input [z, f(x)] +function _distance(v::AbstractVector{T}, s::MOI.Indicator{A}) where {T <: Real} + z = v[1] + # inactive constraint + if A == MOI.ACTIVATE_ON_ONE && isapprox(z, 0) || A == MOI.ACTIVATE_ON_ZERO && isapprox(z, 1) + return zeros(T, 2) + end + return [_distance(z, MOI.ZeroOne()), _distance(v, s.set)] end + +# TODO Complements, requires additional information (bounds for the variables) From 3a786416e73983c3df2c6ed4e0d4e1290848ad8c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Thu, 13 Feb 2020 10:09:39 +0100 Subject: [PATCH 09/28] moved distances to MOI --- src/MathOptInterface.jl | 1 + src/Utilities/Utilities.jl | 1 - .../set_eval.jl => set_distances.jl} | 60 +++++++++---------- 3 files changed, 31 insertions(+), 31 deletions(-) rename src/{Utilities/set_eval.jl => set_distances.jl} (59%) diff --git a/src/MathOptInterface.jl b/src/MathOptInterface.jl index 925bad590c..a1b11fe379 100644 --- a/src/MathOptInterface.jl +++ b/src/MathOptInterface.jl @@ -113,6 +113,7 @@ include("error.jl") include("indextypes.jl") include("functions.jl") include("sets.jl") +include("set_distances.jl") include("attributes.jl") include("constraints.jl") include("modifications.jl") diff --git a/src/Utilities/Utilities.jl b/src/Utilities/Utilities.jl index 1a64e077e5..53c5d3477b 100644 --- a/src/Utilities/Utilities.jl +++ b/src/Utilities/Utilities.jl @@ -32,7 +32,6 @@ include("constraints.jl") include("copy.jl") include("results.jl") include("variables.jl") -include("set_eval.jl") include("model.jl") include("parser.jl") diff --git a/src/Utilities/set_eval.jl b/src/set_distances.jl similarity index 59% rename from src/Utilities/set_eval.jl rename to src/set_distances.jl index 097655d66b..7de0a4d436 100644 --- a/src/Utilities/set_eval.jl +++ b/src/set_distances.jl @@ -8,45 +8,45 @@ function _correct_distance_to_set(d::AbstractVector{T}) where {T <: Real} end function distance_to_set(v, s) - length(v) != MOI.dimension(s) && throw(DimensionMismatch("Mismatch between v and s")) - d = _distance(v, s) + length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between v and s")) + d = unsigned_distance\(v, s) return _correct_distance_to_set(d) end -_distance(v::Real, s::MOI.LessThan) = v - s.upper -_distance(v::Real, s::MOI.GreaterThan) = s.lower - v -_distance(v::Real, s::MOI.EqualTo) = abs(v - s.value) +unsigned_distance\(v::Real, s::LessThan) = v - s.upper +unsigned_distance\(v::Real, s::GreaterThan) = s.lower - v +unsigned_distance\(v::Real, s::EqualTo) = abs(v - s.value) -_distance(v::Real, s::MOI.Interval) = max(s.lower - v, v - s.upper) +unsigned_distance\(v::Real, s::Interval) = max(s.lower - v, v - s.upper) -_distance(v::AbstractVector{<:Real}, ::MOI.Reals) = false * v +unsigned_distance\(v::AbstractVector{<:Real}, ::Reals) = false * v -_distance(v::AbstractVector{<:Real}, ::MOI.Zeros) = abs.(v) +unsigned_distance\(v::AbstractVector{<:Real}, ::Zeros) = abs.(v) -_distance(v::AbstractVector{<:Real}, ::MOI.Nonnegatives) = -v -_distance(v::AbstractVector{<:Real}, ::MOI.Nonpositives) = v +unsigned_distance\(v::AbstractVector{<:Real}, ::Nonnegatives) = -v +unsigned_distance\(v::AbstractVector{<:Real}, ::Nonpositives) = v # Norm cones -function _distance(v::AbstractVector{<:Real}, ::MOI.NormInfinityCone) +function unsigned_distance\(v::AbstractVector{<:Real}, ::NormInfinityCone) t = first(v) xs = v[2:end] return [t - abs(x) for x in xs] end -function _distance(v::AbstractVector{<:Real}, ::MOI.NormOneCone) +function unsigned_distance\(v::AbstractVector{<:Real}, ::NormOneCone) t = v[1] xs = v[2:end] return t - sum(abs, xs) end -function _distance(v::AbstractVector{<:Real}, ::MOI.SecondOrderCone) +function unsigned_distance\(v::AbstractVector{<:Real}, ::SecondOrderCone) t = v[1] xs = v[2:end] return max(-t, dot(xs, xs) - t^2) # avoids sqrt end -function _distance(v::AbstractVector{<:Real}, ::MOI.RotatedSecondOrderCone) +function unsigned_distance\(v::AbstractVector{<:Real}, ::RotatedSecondOrderCone) t = v[1] u = v[2] xs = v[3:end] @@ -56,17 +56,17 @@ function _distance(v::AbstractVector{<:Real}, ::MOI.RotatedSecondOrderCone) ) end -function _distance(v::AbstractVector{<:Real}, s::MOI.GeometricMeanCone) +function unsigned_distance\(v::AbstractVector{<:Real}, s::GeometricMeanCone) t = v[1] xs = v[2:end] - n = MOI.dimension(s) - 1 + n = dimension(s) - 1 return max( -minimum(xs), t^n - prod(xs), ) end -function _distance(v::AbstractVector{<:Real}, ::MOI.ExponentialCone) +function unsigned_distance\(v::AbstractVector{<:Real}, ::ExponentialCone) x = v[1] y = v[2] z = v[3] @@ -76,7 +76,7 @@ function _distance(v::AbstractVector{<:Real}, ::MOI.ExponentialCone) ) end -function _distance(v::AbstractVector{<:Real}, ::MOI.DualExponentialCone) +function unsigned_distance\(v::AbstractVector{<:Real}, ::DualExponentialCone) u = v[1] v = v[2] w = v[3] @@ -86,7 +86,7 @@ function _distance(v::AbstractVector{<:Real}, ::MOI.DualExponentialCone) ) end -function _distance(v::AbstractVector{<:Real}, s::MOI.PowerCone) +function unsigned_distance\(v::AbstractVector{<:Real}, s::PowerCone) x = v[1] y = v[2] z = v[3] @@ -98,7 +98,7 @@ function _distance(v::AbstractVector{<:Real}, s::MOI.PowerCone) ) end -function _distance(v::AbstractVector{<:Real}, s::MOI.DualPowerCone) +function unsigned_distance\(v::AbstractVector{<:Real}, s::DualPowerCone) u = v[1] v = v[2] w = v[3] @@ -111,9 +111,9 @@ function _distance(v::AbstractVector{<:Real}, s::MOI.DualPowerCone) ) end -function _distance(v::AbstractVector{<:Real}, set::MOI.RelativeEntropyCone) +function unsigned_distance\(v::AbstractVector{<:Real}, set::RelativeEntropyCone) all(>=(0), v[2:end]) || return false - n = (MOI.dimension(set)-1) ÷ 2 + n = (dimension(set)-1) ÷ 2 u = v[1] v = v[2:(n+1)] w = v[(n+2):end] @@ -125,14 +125,14 @@ function _distance(v::AbstractVector{<:Real}, set::MOI.RelativeEntropyCone) end -function _distance(v::AbstractVector{<:Real}, s::MOI.NormSpectralCone) +function unsigned_distance\(v::AbstractVector{<:Real}, s::NormSpectralCone) t = v[1] m = reshape(v[2:end], (s.row_dim, s.column_dim)) s1 = LinearAlgebra.svd(m).S[1] return s1 - t end -function _distance(v::AbstractVector{<:Real}, s::MOI.NormNuclearCone) +function unsigned_distance\(v::AbstractVector{<:Real}, s::NormNuclearCone) t = v[1] m = reshape(v[2:end], (s.row_dim, s.column_dim)) s1 = sum(LinearAlgebra.svd(m).S) @@ -141,16 +141,16 @@ end ## Integer sets -function _distance(v::T, ::MOI.ZeroOne) where {T <: Real} +function unsigned_distance\(v::T, ::ZeroOne) where {T <: Real} return min(abs(v - zero(T)), abs(v - one(T))) end -function _distance(v::Real, ::MOI.Integer) +function unsigned_distance\(v::Real, ::Integer) return min(abs(v - floor(v)), abs(v - ceil(v))) end # return the second largest absolute value (=0 if SOS1 respected) -function _distance(v::AbstractVector{T}, s::MOI.SOS1) where {T <: Real} +function unsigned_distance\(v::AbstractVector{T}, s::SOS1) where {T <: Real} first_non_zero = second_non_zero = -1 v1 = zero(T) v2 = zero(T) @@ -182,13 +182,13 @@ end # TODO SOS2 # takes in input [z, f(x)] -function _distance(v::AbstractVector{T}, s::MOI.Indicator{A}) where {T <: Real} +function unsigned_distance\(v::AbstractVector{T}, s::Indicator{A}) where {T <: Real} z = v[1] # inactive constraint - if A == MOI.ACTIVATE_ON_ONE && isapprox(z, 0) || A == MOI.ACTIVATE_ON_ZERO && isapprox(z, 1) + if A == ACTIVATE_ON_ONE && isapprox(z, 0) || A == ACTIVATE_ON_ZERO && isapprox(z, 1) return zeros(T, 2) end - return [_distance(z, MOI.ZeroOne()), _distance(v, s.set)] + return [unsigned_distance\(z, ZeroOne()), unsigned_distance\(v, s.set)] end # TODO Complements, requires additional information (bounds for the variables) From 016fa9da1ed674fc5d7a035c94ac5000c71b293b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Thu, 13 Feb 2020 10:15:27 +0100 Subject: [PATCH 10/28] fix backslash --- src/set_distances.jl | 52 ++++++++++++++++++++++---------------------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/src/set_distances.jl b/src/set_distances.jl index 7de0a4d436..a24d6cebdf 100644 --- a/src/set_distances.jl +++ b/src/set_distances.jl @@ -9,44 +9,44 @@ end function distance_to_set(v, s) length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between v and s")) - d = unsigned_distance\(v, s) + d = unsigned_distance(v, s) return _correct_distance_to_set(d) end -unsigned_distance\(v::Real, s::LessThan) = v - s.upper -unsigned_distance\(v::Real, s::GreaterThan) = s.lower - v -unsigned_distance\(v::Real, s::EqualTo) = abs(v - s.value) +unsigned_distance(v::Real, s::LessThan) = v - s.upper +unsigned_distance(v::Real, s::GreaterThan) = s.lower - v +unsigned_distance(v::Real, s::EqualTo) = abs(v - s.value) -unsigned_distance\(v::Real, s::Interval) = max(s.lower - v, v - s.upper) +unsigned_distance(v::Real, s::Interval) = max(s.lower - v, v - s.upper) -unsigned_distance\(v::AbstractVector{<:Real}, ::Reals) = false * v +unsigned_distance(v::AbstractVector{<:Real}, ::Reals) = false * v -unsigned_distance\(v::AbstractVector{<:Real}, ::Zeros) = abs.(v) +unsigned_distance(v::AbstractVector{<:Real}, ::Zeros) = abs.(v) -unsigned_distance\(v::AbstractVector{<:Real}, ::Nonnegatives) = -v -unsigned_distance\(v::AbstractVector{<:Real}, ::Nonpositives) = v +unsigned_distance(v::AbstractVector{<:Real}, ::Nonnegatives) = -v +unsigned_distance(v::AbstractVector{<:Real}, ::Nonpositives) = v # Norm cones -function unsigned_distance\(v::AbstractVector{<:Real}, ::NormInfinityCone) +function unsigned_distance(v::AbstractVector{<:Real}, ::NormInfinityCone) t = first(v) xs = v[2:end] return [t - abs(x) for x in xs] end -function unsigned_distance\(v::AbstractVector{<:Real}, ::NormOneCone) +function unsigned_distance(v::AbstractVector{<:Real}, ::NormOneCone) t = v[1] xs = v[2:end] return t - sum(abs, xs) end -function unsigned_distance\(v::AbstractVector{<:Real}, ::SecondOrderCone) +function unsigned_distance(v::AbstractVector{<:Real}, ::SecondOrderCone) t = v[1] xs = v[2:end] return max(-t, dot(xs, xs) - t^2) # avoids sqrt end -function unsigned_distance\(v::AbstractVector{<:Real}, ::RotatedSecondOrderCone) +function unsigned_distance(v::AbstractVector{<:Real}, ::RotatedSecondOrderCone) t = v[1] u = v[2] xs = v[3:end] @@ -56,7 +56,7 @@ function unsigned_distance\(v::AbstractVector{<:Real}, ::RotatedSecondOrderCone) ) end -function unsigned_distance\(v::AbstractVector{<:Real}, s::GeometricMeanCone) +function unsigned_distance(v::AbstractVector{<:Real}, s::GeometricMeanCone) t = v[1] xs = v[2:end] n = dimension(s) - 1 @@ -66,7 +66,7 @@ function unsigned_distance\(v::AbstractVector{<:Real}, s::GeometricMeanCone) ) end -function unsigned_distance\(v::AbstractVector{<:Real}, ::ExponentialCone) +function unsigned_distance(v::AbstractVector{<:Real}, ::ExponentialCone) x = v[1] y = v[2] z = v[3] @@ -76,7 +76,7 @@ function unsigned_distance\(v::AbstractVector{<:Real}, ::ExponentialCone) ) end -function unsigned_distance\(v::AbstractVector{<:Real}, ::DualExponentialCone) +function unsigned_distance(v::AbstractVector{<:Real}, ::DualExponentialCone) u = v[1] v = v[2] w = v[3] @@ -86,7 +86,7 @@ function unsigned_distance\(v::AbstractVector{<:Real}, ::DualExponentialCone) ) end -function unsigned_distance\(v::AbstractVector{<:Real}, s::PowerCone) +function unsigned_distance(v::AbstractVector{<:Real}, s::PowerCone) x = v[1] y = v[2] z = v[3] @@ -98,7 +98,7 @@ function unsigned_distance\(v::AbstractVector{<:Real}, s::PowerCone) ) end -function unsigned_distance\(v::AbstractVector{<:Real}, s::DualPowerCone) +function unsigned_distance(v::AbstractVector{<:Real}, s::DualPowerCone) u = v[1] v = v[2] w = v[3] @@ -111,7 +111,7 @@ function unsigned_distance\(v::AbstractVector{<:Real}, s::DualPowerCone) ) end -function unsigned_distance\(v::AbstractVector{<:Real}, set::RelativeEntropyCone) +function unsigned_distance(v::AbstractVector{<:Real}, set::RelativeEntropyCone) all(>=(0), v[2:end]) || return false n = (dimension(set)-1) ÷ 2 u = v[1] @@ -125,14 +125,14 @@ function unsigned_distance\(v::AbstractVector{<:Real}, set::RelativeEntropyCone) end -function unsigned_distance\(v::AbstractVector{<:Real}, s::NormSpectralCone) +function unsigned_distance(v::AbstractVector{<:Real}, s::NormSpectralCone) t = v[1] m = reshape(v[2:end], (s.row_dim, s.column_dim)) s1 = LinearAlgebra.svd(m).S[1] return s1 - t end -function unsigned_distance\(v::AbstractVector{<:Real}, s::NormNuclearCone) +function unsigned_distance(v::AbstractVector{<:Real}, s::NormNuclearCone) t = v[1] m = reshape(v[2:end], (s.row_dim, s.column_dim)) s1 = sum(LinearAlgebra.svd(m).S) @@ -141,16 +141,16 @@ end ## Integer sets -function unsigned_distance\(v::T, ::ZeroOne) where {T <: Real} +function unsigned_distance(v::T, ::ZeroOne) where {T <: Real} return min(abs(v - zero(T)), abs(v - one(T))) end -function unsigned_distance\(v::Real, ::Integer) +function unsigned_distance(v::Real, ::Integer) return min(abs(v - floor(v)), abs(v - ceil(v))) end # return the second largest absolute value (=0 if SOS1 respected) -function unsigned_distance\(v::AbstractVector{T}, s::SOS1) where {T <: Real} +function unsigned_distance(v::AbstractVector{T}, s::SOS1) where {T <: Real} first_non_zero = second_non_zero = -1 v1 = zero(T) v2 = zero(T) @@ -182,13 +182,13 @@ end # TODO SOS2 # takes in input [z, f(x)] -function unsigned_distance\(v::AbstractVector{T}, s::Indicator{A}) where {T <: Real} +function unsigned_distance(v::AbstractVector{T}, s::Indicator{A}) where {T <: Real} 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 [unsigned_distance\(z, ZeroOne()), unsigned_distance\(v, s.set)] + return [unsigned_distance(z, ZeroOne()), unsigned_distance(v, s.set)] end # TODO Complements, requires additional information (bounds for the variables) From 46b6a392f67d29a293698a8bde3ae214e8bb76a8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Thu, 13 Feb 2020 10:19:02 +0100 Subject: [PATCH 11/28] unqualified type param --- src/set_distances.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/set_distances.jl b/src/set_distances.jl index a24d6cebdf..836da75c18 100644 --- a/src/set_distances.jl +++ b/src/set_distances.jl @@ -182,10 +182,10 @@ end # TODO SOS2 # takes in input [z, f(x)] -function unsigned_distance(v::AbstractVector{T}, s::Indicator{A}) where {T <: Real} +function unsigned_distance(v::AbstractVector{T}, s::IndicatorSet{A}) where {A, T <: Real} z = v[1] # inactive constraint - if A == ACTIVATE_ON_ONE && isapprox(z, 0) || A == ACTIVATE_ON_ZERO && isapprox(z, 1) + if A === ACTIVATE_ON_ONE && isapprox(z, 0) || A === ACTIVATE_ON_ZERO && isapprox(z, 1) return zeros(T, 2) end return [unsigned_distance(z, ZeroOne()), unsigned_distance(v, s.set)] From a1d2f6fdc7fa1901e6e865e80ca2c9fcd52d99fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Thu, 13 Feb 2020 11:16:27 +0100 Subject: [PATCH 12/28] documentation, review changes --- src/set_distances.jl | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/src/set_distances.jl b/src/set_distances.jl index 836da75c18..3da84a30b6 100644 --- a/src/set_distances.jl +++ b/src/set_distances.jl @@ -4,11 +4,20 @@ function _correct_distance_to_set(d::Real) end function _correct_distance_to_set(d::AbstractVector{T}) where {T <: Real} - return [ifelse(di < 0, zero(T), di) for di in d] + return [_correct_distance_to_set(di) for di in d] end +""" + 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 `unsigned_distance(v::T, s::S)` with `T` of appropriate type. +""" function distance_to_set(v, s) - length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between v and s")) + length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) d = unsigned_distance(v, s) return _correct_distance_to_set(d) end @@ -19,7 +28,7 @@ unsigned_distance(v::Real, s::EqualTo) = abs(v - s.value) unsigned_distance(v::Real, s::Interval) = max(s.lower - v, v - s.upper) -unsigned_distance(v::AbstractVector{<:Real}, ::Reals) = false * v +unsigned_distance(v::AbstractVector{<:Real}, ::Reals) = zeros(length(v)) unsigned_distance(v::AbstractVector{<:Real}, ::Zeros) = abs.(v) From 90d36314437742b3ce6d761c789b6a1f41f26976 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Thu, 13 Feb 2020 12:30:53 +0100 Subject: [PATCH 13/28] vectorized versions --- src/set_distances.jl | 77 ++++++++++++-------------------------------- 1 file changed, 20 insertions(+), 57 deletions(-) diff --git a/src/set_distances.jl b/src/set_distances.jl index 3da84a30b6..0ea817e224 100644 --- a/src/set_distances.jl +++ b/src/set_distances.jl @@ -52,47 +52,35 @@ end function unsigned_distance(v::AbstractVector{<:Real}, ::SecondOrderCone) t = v[1] xs = v[2:end] - return max(-t, dot(xs, xs) - t^2) # avoids sqrt + return [-t, dot(xs, xs) - t^2] # avoids sqrt end function unsigned_distance(v::AbstractVector{<:Real}, ::RotatedSecondOrderCone) t = v[1] u = v[2] xs = v[3:end] - return max( - -t, -u, - dot(xs, xs) - 2 * t * u - ) + return [-t, -u, dot(xs, xs) - 2 * t * u] end function unsigned_distance(v::AbstractVector{<:Real}, s::GeometricMeanCone) t = v[1] xs = v[2:end] n = dimension(s) - 1 - return max( - -minimum(xs), - t^n - prod(xs), - ) + return [-minimum(xs), t^n - prod(xs)] end function unsigned_distance(v::AbstractVector{<:Real}, ::ExponentialCone) x = v[1] y = v[2] z = v[3] - return max( - y, - y * exp(x/y) - z, - ) + return [y, y * exp(x/y) - z] end function unsigned_distance(v::AbstractVector{<:Real}, ::DualExponentialCone) u = v[1] v = v[2] w = v[3] - return max( - -u, - -u*exp(v/u) - ℯ * w - ) + return [-u, -u*exp(v/u) - ℯ * w] end function unsigned_distance(v::AbstractVector{<:Real}, s::PowerCone) @@ -100,11 +88,7 @@ function unsigned_distance(v::AbstractVector{<:Real}, s::PowerCone) y = v[2] z = v[3] e = s.exponent - return max( - x, - y, - abs(z) - x^e * y^(1-e) - ) + return [-x, -y, abs(z) - x^e * y^(1-e)] end function unsigned_distance(v::AbstractVector{<:Real}, s::DualPowerCone) @@ -113,11 +97,7 @@ function unsigned_distance(v::AbstractVector{<:Real}, s::DualPowerCone) w = v[3] e = s.exponent ce = 1-e - return max( - u, - v, - abs(w) - (u/e)^e * (v/ce)^ce - ) + return [u, v, abs(w) - (u/e)^e * (v/ce)^ce] end function unsigned_distance(v::AbstractVector{<:Real}, set::RelativeEntropyCone) @@ -127,10 +107,7 @@ function unsigned_distance(v::AbstractVector{<:Real}, set::RelativeEntropyCone) v = v[2:(n+1)] w = v[(n+2):end] s = sum(w[i] * log(w[i]/v[i]) for i in eachindex(w)) - return max( - -minimum(v[2:end]), - s - u, - ) + return [-minimum(v[2:end]), s - u] end @@ -158,34 +135,20 @@ function unsigned_distance(v::Real, ::Integer) return min(abs(v - floor(v)), abs(v - ceil(v))) end -# return the second largest absolute value (=0 if SOS1 respected) -function unsigned_distance(v::AbstractVector{T}, s::SOS1) where {T <: Real} - first_non_zero = second_non_zero = -1 - v1 = zero(T) - v2 = zero(T) - for (i, x) in enumerate(v) - if !isapprox(x, 0) - if first_non_zero < 0 - first_non_zero = i - v1 = x - elseif second_non_zero < 0 - second_non_zero = i - v2 = x - else # replace one of the two - xa = abs(x) - if xa > v1 - second_non_zero = first_non_zero - v2 = v1 - first_non_zero = i - v1 = xa - elseif v[second_non_zero] <= x < v[first_non_zero] - second_non_zero = i - v2 = xa - end - end +# return the element-wise distance to zero, with the greatest element to 0 +function unsigned_distance(v::AbstractVector{T}, ::SOS1) where {T <: Real} + d = unsigned_distance(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 - return v2 end # TODO SOS2 From ae2c2a84343e28f421b508acb4bcfcf541f407e9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Thu, 13 Feb 2020 22:16:48 +0100 Subject: [PATCH 14/28] add boilerplate to all methods --- src/set_distances.jl | 107 +++++++++++++++++++++++++------------------ 1 file changed, 62 insertions(+), 45 deletions(-) diff --git a/src/set_distances.jl b/src/set_distances.jl index 0ea817e224..9fc2ce3127 100644 --- a/src/set_distances.jl +++ b/src/set_distances.jl @@ -14,130 +14,147 @@ 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 `unsigned_distance(v::T, s::S)` with `T` of appropriate type. +Each set `S` implements `distance_to_set(v::T, s::S)` with `T` of appropriate type. """ -function distance_to_set(v, s) - length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) - d = unsigned_distance(v, s) - return _correct_distance_to_set(d) -end +function distance_to_set end -unsigned_distance(v::Real, s::LessThan) = v - s.upper -unsigned_distance(v::Real, s::GreaterThan) = s.lower - v -unsigned_distance(v::Real, s::EqualTo) = abs(v - s.value) +distance_to_set(v::Real, s::LessThan) = _correct_distance_to_set(v - s.upper) +distance_to_set(v::Real, s::GreaterThan) = _correct_distance_to_set(s.lower - v) +distance_to_set(v::Real, s::EqualTo) = abs(v - s.value) -unsigned_distance(v::Real, s::Interval) = max(s.lower - v, v - s.upper) +distance_to_set(v::T, s::Interval) where {T <: Real} = max(s.lower - v, v - s.upper, zero(T)) -unsigned_distance(v::AbstractVector{<:Real}, ::Reals) = zeros(length(v)) +function distance_to_set(v::AbstractVector{<:Real}, s::Reals) + length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) + return zeros(length(v)) +end -unsigned_distance(v::AbstractVector{<:Real}, ::Zeros) = abs.(v) +function distance_to_set(v::AbstractVector{<:Real}, s::Zeros) + length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) + return abs.(v) +end -unsigned_distance(v::AbstractVector{<:Real}, ::Nonnegatives) = -v -unsigned_distance(v::AbstractVector{<:Real}, ::Nonpositives) = v +function distance_to_set(v::AbstractVector{<:Real}, s::Nonnegatives) + length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) + return _correct_distance_to_set(-v) +end +function distance_to_set(v::AbstractVector{<:Real}, s::Nonpositives) + length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) + return _correct_distance_to_set(v) +end # Norm cones -function unsigned_distance(v::AbstractVector{<:Real}, ::NormInfinityCone) +function distance_to_set(v::AbstractVector{<:Real}, s::NormInfinityCone) + length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) t = first(v) xs = v[2:end] - return [t - abs(x) for x in xs] + return [max(t - abs(x), zero(t)) for x in xs] end -function unsigned_distance(v::AbstractVector{<:Real}, ::NormOneCone) +function distance_to_set(v::AbstractVector{<:Real}, s::NormOneCone) + length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) t = v[1] xs = v[2:end] - return t - sum(abs, xs) + return _correct_distance_to_set(t - sum(abs, xs)) end -function unsigned_distance(v::AbstractVector{<:Real}, ::SecondOrderCone) +function distance_to_set(v::AbstractVector{<:Real}, s::SecondOrderCone) + length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) t = v[1] xs = v[2:end] - return [-t, dot(xs, xs) - t^2] # avoids sqrt + return _correct_distance_to_set([-t, dot(xs, xs) - t^2]) # avoids sqrt end -function unsigned_distance(v::AbstractVector{<:Real}, ::RotatedSecondOrderCone) +function distance_to_set(v::AbstractVector{<:Real}, s::RotatedSecondOrderCone) + length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) t = v[1] u = v[2] xs = v[3:end] - return [-t, -u, dot(xs, xs) - 2 * t * u] + return _correct_distance_to_set([-t, -u, dot(xs, xs) - 2 * t * u]) end -function unsigned_distance(v::AbstractVector{<:Real}, s::GeometricMeanCone) +function distance_to_set(v::AbstractVector{<:Real}, s::GeometricMeanCone) + length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) t = v[1] xs = v[2:end] n = dimension(s) - 1 - return [-minimum(xs), t^n - prod(xs)] + return _correct_distance_to_set([-minimum(xs), t^n - prod(xs)]) end -function unsigned_distance(v::AbstractVector{<:Real}, ::ExponentialCone) +function distance_to_set(v::AbstractVector{<:Real}, s::ExponentialCone) + length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) x = v[1] y = v[2] z = v[3] - return [y, y * exp(x/y) - z] + return _correct_distance_to_set([y, y * exp(x/y) - z]) end -function unsigned_distance(v::AbstractVector{<:Real}, ::DualExponentialCone) +function distance_to_set(v::AbstractVector{<:Real}, s::DualExponentialCone) + length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) u = v[1] v = v[2] w = v[3] - return [-u, -u*exp(v/u) - ℯ * w] + return _correct_distance_to_set([-u, -u*exp(v/u) - ℯ * w]) end -function unsigned_distance(v::AbstractVector{<:Real}, s::PowerCone) +function distance_to_set(v::AbstractVector{<:Real}, s::PowerCone) + length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) x = v[1] y = v[2] z = v[3] e = s.exponent - return [-x, -y, abs(z) - x^e * y^(1-e)] + return _correct_distance_to_set([-x, -y, abs(z) - x^e * y^(1-e)]) end -function unsigned_distance(v::AbstractVector{<:Real}, s::DualPowerCone) +function distance_to_set(v::AbstractVector{<:Real}, s::DualPowerCone) + length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) u = v[1] v = v[2] w = v[3] e = s.exponent ce = 1-e - return [u, v, abs(w) - (u/e)^e * (v/ce)^ce] + return _correct_distance_to_set([u, v, abs(w) - (u/e)^e * (v/ce)^ce]) end -function unsigned_distance(v::AbstractVector{<:Real}, set::RelativeEntropyCone) +function distance_to_set(v::AbstractVector{<:Real}, set::RelativeEntropyCone) + length(v) != dimension(set) && throw(DimensionMismatch("Mismatch between value and set")) all(>=(0), v[2:end]) || return false 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)) - return [-minimum(v[2:end]), s - u] + return _correct_distance_to_set([-minimum(v[2:end]), s - u]) end - -function unsigned_distance(v::AbstractVector{<:Real}, s::NormSpectralCone) +function distance_to_set(v::AbstractVector{<:Real}, s::NormSpectralCone) t = v[1] m = reshape(v[2:end], (s.row_dim, s.column_dim)) s1 = LinearAlgebra.svd(m).S[1] - return s1 - t + return _correct_distance_to_set(s1 - t) end -function unsigned_distance(v::AbstractVector{<:Real}, s::NormNuclearCone) +function distance_to_set(v::AbstractVector{<:Real}, s::NormNuclearCone) t = v[1] m = reshape(v[2:end], (s.row_dim, s.column_dim)) s1 = sum(LinearAlgebra.svd(m).S) - return s1 - t + return _correct_distance_to_set(s1 - t) end ## Integer sets -function unsigned_distance(v::T, ::ZeroOne) where {T <: Real} +function distance_to_set(v::T, ::ZeroOne) where {T <: Real} return min(abs(v - zero(T)), abs(v - one(T))) end -function unsigned_distance(v::Real, ::Integer) +function distance_to_set(v::Real, ::Integer) return min(abs(v - floor(v)), abs(v - ceil(v))) end # return the element-wise distance to zero, with the greatest element to 0 -function unsigned_distance(v::AbstractVector{T}, ::SOS1) where {T <: Real} - d = unsigned_distance(v, Zeros(length(v))) +function distance_to_set(v::AbstractVector{T}, ::SOS1) where {T <: Real} + d = distance_to_set(v, Zeros(length(v))) m = maximum(d) if m ≈ zero(T) return d @@ -154,13 +171,13 @@ end # TODO SOS2 # takes in input [z, f(x)] -function unsigned_distance(v::AbstractVector{T}, s::IndicatorSet{A}) where {A, T <: Real} +function distance_to_set(v::AbstractVector{T}, s::IndicatorSet{A}) where {A, T <: Real} 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 [unsigned_distance(z, ZeroOne()), unsigned_distance(v, s.set)] + return [distance_to_set(z, ZeroOne()), distance_to_set(v, s.set)] end # TODO Complements, requires additional information (bounds for the variables) From b3678b587791287c44e0683db5caf3dc372a3a20 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Thu, 20 Feb 2020 10:24:12 +0100 Subject: [PATCH 15/28] check all max locally --- src/set_distances.jl | 97 ++++++++++++++++++++++++++------------------ 1 file changed, 58 insertions(+), 39 deletions(-) diff --git a/src/set_distances.jl b/src/set_distances.jl index 9fc2ce3127..474e3fdcd8 100644 --- a/src/set_distances.jl +++ b/src/set_distances.jl @@ -1,12 +1,3 @@ - -function _correct_distance_to_set(d::Real) - return ifelse(d < 0, zero(d), d) -end - -function _correct_distance_to_set(d::AbstractVector{T}) where {T <: Real} - return [_correct_distance_to_set(di) for di in d] -end - """ distance_to_set(v, s) @@ -18,128 +9,154 @@ Each set `S` implements `distance_to_set(v::T, s::S)` with `T` of appropriate ty """ function distance_to_set end -distance_to_set(v::Real, s::LessThan) = _correct_distance_to_set(v - s.upper) -distance_to_set(v::Real, s::GreaterThan) = _correct_distance_to_set(s.lower - v) +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) distance_to_set(v::T, s::Interval) where {T <: Real} = max(s.lower - v, v - s.upper, zero(T)) -function distance_to_set(v::AbstractVector{<:Real}, s::Reals) +function _check_dimension(v::AbstractVector, s) length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) + return nothing +end + +function distance_to_set(v::AbstractVector{<:Real}, s::Reals) + _check_dimension(v, s) return zeros(length(v)) end function distance_to_set(v::AbstractVector{<:Real}, s::Zeros) - length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) + _check_dimension(v, s) return abs.(v) end function distance_to_set(v::AbstractVector{<:Real}, s::Nonnegatives) - length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) - return _correct_distance_to_set(-v) + _check_dimension(v, s) + return [ifelse(vi < 0, -vi, zero(vi)) for vi in v] end function distance_to_set(v::AbstractVector{<:Real}, s::Nonpositives) - length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) - return _correct_distance_to_set(v) + _check_dimension(v, s) + return [ifelse(vi > 0, vi, zero(vi)) for vi in v] end # Norm cones function distance_to_set(v::AbstractVector{<:Real}, s::NormInfinityCone) - length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) + _check_dimension(v, s) t = first(v) xs = v[2:end] return [max(t - abs(x), zero(t)) for x in xs] end function distance_to_set(v::AbstractVector{<:Real}, s::NormOneCone) - length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) + _check_dimension(v, s) t = v[1] xs = v[2:end] - return _correct_distance_to_set(t - sum(abs, xs)) + result = t - sum(abs, xs) + return max(result, zero(result)) end function distance_to_set(v::AbstractVector{<:Real}, s::SecondOrderCone) - length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) + _check_dimension(v, s) t = v[1] xs = v[2:end] - return _correct_distance_to_set([-t, dot(xs, xs) - t^2]) # avoids sqrt + result = LinearAlgebra.norm(xs) - t + return max(result, zero(result)) # avoids sqrt end function distance_to_set(v::AbstractVector{<:Real}, s::RotatedSecondOrderCone) - length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) + _check_dimension(v, s) t = v[1] u = v[2] xs = v[3:end] - return _correct_distance_to_set([-t, -u, dot(xs, xs) - 2 * t * u]) + return [ + max(-t, zero(t)), + max(-u, zero(u)), + max(dot(xs,xs) - 2 * t * u), + ] end function distance_to_set(v::AbstractVector{<:Real}, s::GeometricMeanCone) - length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) + _check_dimension(v, s) t = v[1] xs = v[2:end] n = dimension(s) - 1 - return _correct_distance_to_set([-minimum(xs), t^n - prod(xs)]) + result = t^n - prod(xs) + return push!( + max.(xs, zero(result)), # x >= 0 + max(result, zero(result)), + ) end function distance_to_set(v::AbstractVector{<:Real}, s::ExponentialCone) - length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) + _check_dimension(v, s) x = v[1] y = v[2] z = v[3] - return _correct_distance_to_set([y, y * exp(x/y) - z]) + result = y * exp(x/y) - z + return [max(y, zero(result)), max(result, zero(result))] end function distance_to_set(v::AbstractVector{<:Real}, s::DualExponentialCone) - length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) + _check_dimension(v, s) u = v[1] v = v[2] w = v[3] - return _correct_distance_to_set([-u, -u*exp(v/u) - ℯ * w]) + result = -u*exp(v/u) - ℯ * w + return [max(-u, zero(result)), max(result, zero(result))] end function distance_to_set(v::AbstractVector{<:Real}, s::PowerCone) - length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) + _check_dimension(v, s) x = v[1] y = v[2] z = v[3] e = s.exponent - return _correct_distance_to_set([-x, -y, abs(z) - x^e * y^(1-e)]) + result = abs(z) - x^e * y^(1-e) + return [max(-x, zero(result)), max(-y, zero(result)), max(result, zero(result))] end function distance_to_set(v::AbstractVector{<:Real}, s::DualPowerCone) - length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) + _check_dimension(v, s) u = v[1] v = v[2] w = v[3] e = s.exponent ce = 1-e - return _correct_distance_to_set([u, v, abs(w) - (u/e)^e * (v/ce)^ce]) + result = abs(w) - (u/e)^e * (v/ce)^ce + return [max(-u, zero(result)), max(-v, zero(result)), max(result, zero(result))] end function distance_to_set(v::AbstractVector{<:Real}, set::RelativeEntropyCone) - length(v) != dimension(set) && throw(DimensionMismatch("Mismatch between value and set")) - all(>=(0), v[2:end]) || return false + _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)) - return _correct_distance_to_set([-minimum(v[2:end]), s - u]) + result = s - u + return push!( + max.(v[2:end], zero(result)), + max(result, zero(result)), + ) end 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] - return _correct_distance_to_set(s1 - t) + result = s1 - t + return max(result, zero(result)) end 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) - return _correct_distance_to_set(s1 - t) + result = s1 - t + return max(result, zero(result)) end ## Integer sets @@ -154,6 +171,7 @@ 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) @@ -172,6 +190,7 @@ end # 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) From 73aaa3b1ead143a71d04ed0f22e19eb8ca28c1b7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Thu, 20 Feb 2020 10:36:55 +0100 Subject: [PATCH 16/28] set distance in sets --- src/MathOptInterface.jl | 2 - src/set_distances.jl | 202 ---------------------------------------- src/sets.jl | 189 +++++++++++++++++++++++++++++++++++++ 3 files changed, 189 insertions(+), 204 deletions(-) delete mode 100644 src/set_distances.jl diff --git a/src/MathOptInterface.jl b/src/MathOptInterface.jl index a1b11fe379..6d64399309 100644 --- a/src/MathOptInterface.jl +++ b/src/MathOptInterface.jl @@ -113,8 +113,6 @@ include("error.jl") include("indextypes.jl") include("functions.jl") include("sets.jl") -include("set_distances.jl") -include("attributes.jl") include("constraints.jl") include("modifications.jl") include("variables.jl") diff --git a/src/set_distances.jl b/src/set_distances.jl deleted file mode 100644 index 474e3fdcd8..0000000000 --- a/src/set_distances.jl +++ /dev/null @@ -1,202 +0,0 @@ -""" - 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 - -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) - -distance_to_set(v::T, s::Interval) where {T <: Real} = max(s.lower - v, v - s.upper, zero(T)) - -function _check_dimension(v::AbstractVector, s) - length(v) != dimension(s) && throw(DimensionMismatch("Mismatch between value and set")) - return nothing -end - -function distance_to_set(v::AbstractVector{<:Real}, s::Reals) - _check_dimension(v, s) - return zeros(length(v)) -end - -function distance_to_set(v::AbstractVector{<:Real}, s::Zeros) - _check_dimension(v, s) - return abs.(v) -end - -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 -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 - -# Norm cones - -function distance_to_set(v::AbstractVector{<:Real}, s::NormInfinityCone) - _check_dimension(v, s) - t = first(v) - xs = v[2:end] - return [max(t - abs(x), zero(t)) for x in xs] -end - -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 - -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 - -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 - -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 - -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 - -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 - -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 - -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 - -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 - -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 - -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 - -## Integer sets - -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 - -# 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 - -# TODO SOS2 - -# 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, s.set)] -end - -# TODO Complements, requires additional information (bounds for the variables) diff --git a/src/sets.jl b/src/sets.jl index 4bbfcc5ee8..38742b8163 100644 --- a/src/sets.jl +++ b/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 From 3cab392d2f512ce159ba8d12de1f12fab1a340ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Fri, 21 Feb 2020 09:15:27 +0100 Subject: [PATCH 17/28] removed file --- src/MathOptInterface.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/MathOptInterface.jl b/src/MathOptInterface.jl index 6d64399309..925bad590c 100644 --- a/src/MathOptInterface.jl +++ b/src/MathOptInterface.jl @@ -113,6 +113,7 @@ include("error.jl") include("indextypes.jl") include("functions.jl") include("sets.jl") +include("attributes.jl") include("constraints.jl") include("modifications.jl") include("variables.jl") From 83a6131fad0b5a15b9df37f3b2d0174d1a720deb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Fri, 6 Mar 2020 15:24:32 +0100 Subject: [PATCH 18/28] norm on all multi-results --- src/sets.jl | 39 +++++++++++++++++++++++++-------------- 1 file changed, 25 insertions(+), 14 deletions(-) diff --git a/src/sets.jl b/src/sets.jl index 38742b8163..f451dfee85 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -128,9 +128,9 @@ end dual_set(s::Reals) = Zeros(dimension(s)) dual_set_type(::Type{Reals}) = Zeros -function distance_to_set(v::AbstractVector{<:Real}, s::Reals) +function distance_to_set(v::AbstractVector{T}, s::Reals) where {T <: Real} _check_dimension(v, s) - return zeros(length(v)) + return zero(T) end """ @@ -147,7 +147,7 @@ dual_set_type(::Type{Zeros}) = Reals function distance_to_set(v::AbstractVector{<:Real}, s::Zeros) _check_dimension(v, s) - return abs.(v) + return LinearAlgebra.norm2(v) end """ @@ -164,7 +164,7 @@ 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] + return LinearAlgebra.norm2(ifelse(vi < 0, -vi, zero(vi)) for vi in v) end """ @@ -181,7 +181,7 @@ 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] + return LinearAlgebra.norm2(ifelse(vi > 0, vi, zero(vi)) for vi in v) end """ @@ -348,10 +348,11 @@ function distance_to_set(v::AbstractVector{<:Real}, s::GeometricMeanCone) xs = v[2:end] n = dimension(s) - 1 result = t^n - prod(xs) - return push!( + return + LinearAlgebra.norm2(push!( max.(xs, zero(result)), # x >= 0 max(result, zero(result)), - ) + )) end """ @@ -370,7 +371,9 @@ function distance_to_set(v::AbstractVector{<:Real}, s::ExponentialCone) y = v[2] z = v[3] result = y * exp(x/y) - z - return [max(y, zero(result)), max(result, zero(result))] + return LinearAlgebra.norm2( + (max(y, zero(result)), max(result, zero(result))) + ) end """ @@ -389,7 +392,9 @@ function distance_to_set(v::AbstractVector{<:Real}, s::DualExponentialCone) v = v[2] w = v[3] result = -u*exp(v/u) - ℯ * w - return [max(-u, zero(result)), max(result, zero(result))] + return LinearAlgebra.norm2( + (max(-u, zero(result)), max(result, zero(result))) + ) end """ @@ -411,7 +416,9 @@ function distance_to_set(v::AbstractVector{<:Real}, s::PowerCone) 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))] + return LinearAlgebra.norm2( + (max(-x, zero(result)), max(-y, zero(result)), max(result, zero(result))) + ) end """ @@ -434,7 +441,9 @@ function distance_to_set(v::AbstractVector{<:Real}, s::DualPowerCone) 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))] + return LinearAlgebra.norm2( + (max(-u, zero(result)), max(-v, zero(result)), max(result, zero(result))) + ) end dimension(s::Union{ExponentialCone, DualExponentialCone, PowerCone, DualPowerCone}) = 3 @@ -460,10 +469,10 @@ function distance_to_set(v::AbstractVector{<:Real}, set::RelativeEntropyCone) w = v[(n+2):end] s = sum(w[i] * log(w[i]/v[i]) for i in eachindex(w)) result = s - u - return push!( + return LinearAlgebra.norm2(push!( max.(v[2:end], zero(result)), max(result, zero(result)), - ) + )) end """ @@ -930,7 +939,9 @@ function distance_to_set(v::AbstractVector{T}, s::IndicatorSet{A}) where {A, T < 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)] + return LinearAlgebra.norm2( + (distance_to_set(z, ZeroOne()), distance_to_set(v[2], s.set)) + ) end function Base.copy(set::IndicatorSet{A,S}) where {A,S} From c9690689b7c1e66ab116137d51e5b8c8e9fc4d5c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Fri, 6 Mar 2020 18:40:17 +0100 Subject: [PATCH 19/28] first tests --- src/sets.jl | 37 +++++++++++++++-------- test/runtests.jl | 1 + test/set_distances.jl | 69 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 94 insertions(+), 13 deletions(-) create mode 100644 test/set_distances.jl diff --git a/src/sets.jl b/src/sets.jl index f451dfee85..dc8102fe9f 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -1,6 +1,8 @@ # Sets # Note: When adding a new set, also add it to Utilities.Model. +import LinearAlgebra +using LinearAlgebra: dot """ AbstractSet @@ -269,6 +271,14 @@ end dual_set(s::NormInfinityCone) = NormOneCone(dimension(s)) dual_set_type(::Type{NormInfinityCone}) = NormOneCone +function distance_to_set(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) @@ -285,7 +295,7 @@ 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) + result = sum(abs, xs) - t return max(result, zero(result)) end @@ -305,7 +315,7 @@ 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 + result = LinearAlgebra.norm2(xs) - t return max(result, zero(result)) # avoids sqrt end @@ -326,11 +336,9 @@ function distance_to_set(v::AbstractVector{<:Real}, s::RotatedSecondOrderCone) 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), - ] + return LinearAlgebra.norm2( + (max(-t, zero(t)), max(-u, zero(u)), max(dot(xs,xs) - 2 * t * u)) + ) end """ @@ -347,12 +355,15 @@ function distance_to_set(v::AbstractVector{<:Real}, s::GeometricMeanCone) t = v[1] xs = v[2:end] n = dimension(s) - 1 - result = t^n - prod(xs) - return - LinearAlgebra.norm2(push!( - max.(xs, zero(result)), # x >= 0 - max(result, zero(result)), - )) + 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 """ diff --git a/test/runtests.jl b/test/runtests.jl index 6c530b84ad..6939134e42 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,6 +21,7 @@ include("dummy.jl") include("errors.jl") include("functions.jl") include("sets.jl") + include("set_distances.jl") include("attributes.jl") include("instantiate.jl") end diff --git a/test/set_distances.jl b/test/set_distances.jl new file mode 100644 index 0000000000..667806e51d --- /dev/null +++ b/test/set_distances.jl @@ -0,0 +1,69 @@ +using Test + +using MathOptInterface +const MOI = MathOptInterface +import LinearAlgebra + +@testset "Set distances" begin + @testset "$n-dimensional orthants" for n in 1:3:15 + v = rand(n) + @test MOI.distance_to_set(v, MOI.Reals(n)) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(v, MOI.Zeros(n)) > 0 + @test MOI.distance_to_set(v, MOI.Nonnegatives(n)) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(-v, MOI.Nonpositives(n)) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(v, MOI.Nonpositives(n)) ≈ MOI.distance_to_set(-v, MOI.Nonnegatives(n)) > 0 + end + + @testset "Scalar comparisons" begin + values = rand(10) + for v in values + @test MOI.distance_to_set(v, MOI.EqualTo(v)) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(-v, MOI.EqualTo(v)) ≈ MOI.distance_to_set(v, MOI.EqualTo(-v)) ≈ 2v + @test MOI.distance_to_set(v, MOI.LessThan(v)) ≈ MOI.distance_to_set(v, MOI.LessThan(v+1)) ≈ 0 + @test MOI.distance_to_set(v, MOI.LessThan(0)) ≈ MOI.distance_to_set(-v, MOI.GreaterThan(0)) ≈ v + @test MOI.distance_to_set(v, MOI.GreaterThan(v)) ≈ MOI.distance_to_set(v+1, MOI.GreaterThan(v+1)) ≈ 0 + @test MOI.distance_to_set(v, MOI.Interval(v,v)) ≈ MOI.distance_to_set(v, MOI.Interval(-v,v)) ≈ 0 + @test MOI.distance_to_set(v, MOI.Interval(-v, 0.0)) ≈ MOI.distance_to_set(-v, MOI.Interval(0.0, v)) ≈ v + end + end + @testset "$n-dimensional norm cones" for n in 2:5:15 + x = rand(n) + tsum = sum(x) + vsum = vcat(tsum, x) + @test MOI.distance_to_set(vsum, MOI.NormOneCone(n+1)) ≈ MOI.distance_to_set(vsum, MOI.NormInfinityCone(n+1)) ≈ 0 + tmax = maximum(x) + vmax = vcat(tmax, x) + @test MOI.distance_to_set(vmax, MOI.NormOneCone(n+1)) > 0 + @test MOI.distance_to_set(vmax, MOI.NormInfinityCone(n+1)) ≈ 0 atol=eps(Float64) + tmin = 0 + vmin = vcat(tmin, x) + @test MOI.distance_to_set(vmin, MOI.NormInfinityCone(n+1)) ≈ tmax + @test MOI.distance_to_set(vmin, MOI.NormOneCone(n+1)) ≈ tsum + + tvalid = sqrt(n) # upper bound on the norm2 + vok_soc = vcat(tvalid, x) + @test MOI.distance_to_set(vok_soc, MOI.SecondOrderCone(n+1) ) ≈ 0 atol=eps(Float64) + vko_soc = vcat(-2, x) + @test MOI.distance_to_set(vko_soc, MOI.SecondOrderCone(n+1) ) ≈ 2 + LinearAlgebra.norm2(x) + + vko_soc = vcat(-2, x) + @test MOI.distance_to_set(vko_soc, MOI.SecondOrderCone(n+1) ) ≈ 2 + LinearAlgebra.norm2(x) + + t_ko_rot = u_ko_rot = LinearAlgebra.norm2(x) / 2 + vko_roc = vcat(t_ko_rot, u_ko_rot, x) + @test MOI.distance_to_set(vko_roc, MOI.RotatedSecondOrderCone(n+2)) ≈ LinearAlgebra.dot(x,x) / 2 + vok_roc = vcat(t_ko_rot * 2, u_ko_rot, x) + @test MOI.distance_to_set(vok_roc, MOI.RotatedSecondOrderCone(n+2)) ≈ 0 atol=eps(Float64) + end + + @testset "Other vector cones" for n in 2:5:15 + x = rand(n) + t = 0.5 * prod(x)^(inv(n)) + vok = vcat(t, x) + @test MOI.distance_to_set(vok, MOI.GeometricMeanCone(n+1)) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(vcat(t / 2, x), MOI.GeometricMeanCone(n+1)) ≈ 0 atol=eps(Float64) + # negative x always means positive distance + @test MOI.distance_to_set(vcat(t / 2, vcat(x, -1)), MOI.GeometricMeanCone(n+2)) > 0 + @test MOI.distance_to_set(vcat(t / 2, -x), MOI.GeometricMeanCone(n+1)) > 0 + end +end From 891426ff60ec4d66c625f1adf11f9a286a5347a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Fri, 6 Mar 2020 19:09:55 +0100 Subject: [PATCH 20/28] exponential cones --- src/sets.jl | 12 ++++++------ test/set_distances.jl | 24 +++++++++++++++++++++--- 2 files changed, 27 insertions(+), 9 deletions(-) diff --git a/src/sets.jl b/src/sets.jl index dc8102fe9f..fc48f43958 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -383,7 +383,7 @@ function distance_to_set(v::AbstractVector{<:Real}, s::ExponentialCone) z = v[3] result = y * exp(x/y) - z return LinearAlgebra.norm2( - (max(y, zero(result)), max(result, zero(result))) + (max(-y, zero(result)), max(result, zero(result))) ) end @@ -397,11 +397,11 @@ 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] +function distance_to_set(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))) diff --git a/test/set_distances.jl b/test/set_distances.jl index 667806e51d..9216824e97 100644 --- a/test/set_distances.jl +++ b/test/set_distances.jl @@ -53,10 +53,10 @@ import LinearAlgebra vko_roc = vcat(t_ko_rot, u_ko_rot, x) @test MOI.distance_to_set(vko_roc, MOI.RotatedSecondOrderCone(n+2)) ≈ LinearAlgebra.dot(x,x) / 2 vok_roc = vcat(t_ko_rot * 2, u_ko_rot, x) - @test MOI.distance_to_set(vok_roc, MOI.RotatedSecondOrderCone(n+2)) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(vok_roc, MOI.RotatedSecondOrderCone(n+2)) ≈ 0 atol=5eps(Float64) end - @testset "Other vector cones" for n in 2:5:15 + @testset "Geometric Mean cone dimension $n" for n in 2:5:15 x = rand(n) t = 0.5 * prod(x)^(inv(n)) vok = vcat(t, x) @@ -64,6 +64,24 @@ import LinearAlgebra @test MOI.distance_to_set(vcat(t / 2, x), MOI.GeometricMeanCone(n+1)) ≈ 0 atol=eps(Float64) # negative x always means positive distance @test MOI.distance_to_set(vcat(t / 2, vcat(x, -1)), MOI.GeometricMeanCone(n+2)) > 0 - @test MOI.distance_to_set(vcat(t / 2, -x), MOI.GeometricMeanCone(n+1)) > 0 + @test MOI.distance_to_set(vcat(t / 2, -x), MOI.GeometricMeanCone(n+1)) > 0 + end + @testset "Exponential cones" for _ in 1:30 + (x, y, z) = rand(3) + y += 1 # ensure y > 0 + if y * exp(x/y) <= z + @test MOI.distance_to_set([x, y, z], MOI.ExponentialCone()) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set([x, -1, z], MOI.ExponentialCone()) ≈ 1 atol=eps(Float64) + else + @test MOI.distance_to_set([x, y, z], MOI.ExponentialCone()) ≈ y * exp(x/y) - z + end + (u, v, w) = randn(3) + if u != 0.0 # just in case not to blow up + if u*exp(v/u) < ℯ * w && u < 0 + @test MOI.distance_to_set([u, v, w], MOI.DualExponentialCone()) ≈ 0 atol=eps(Float64) + elseif u < 0 + @test MOI.distance_to_set([u, v, w], MOI.DualExponentialCone()) ≈ u*exp(v/u) - ℯ * w + end + end end end From 5cde8c80b26190e61ccad39cdc9e782c78399c6e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Fri, 6 Mar 2020 19:30:54 +0100 Subject: [PATCH 21/28] power cone --- src/sets.jl | 8 +++++++- test/set_distances.jl | 20 +++++++++++++++++++- 2 files changed, 26 insertions(+), 2 deletions(-) diff --git a/src/sets.jl b/src/sets.jl index fc48f43958..1aa6b5709b 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -426,9 +426,15 @@ function distance_to_set(v::AbstractVector{<:Real}, s::PowerCone) 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(result)), max(-y, zero(result))) + ) + end result = abs(z) - x^e * y^(1-e) return LinearAlgebra.norm2( - (max(-x, zero(result)), max(-y, zero(result)), max(result, zero(result))) + max(result, zero(result)) ) end diff --git a/test/set_distances.jl b/test/set_distances.jl index 9216824e97..c1f41dd1b9 100644 --- a/test/set_distances.jl +++ b/test/set_distances.jl @@ -66,7 +66,8 @@ import LinearAlgebra @test MOI.distance_to_set(vcat(t / 2, vcat(x, -1)), MOI.GeometricMeanCone(n+2)) > 0 @test MOI.distance_to_set(vcat(t / 2, -x), MOI.GeometricMeanCone(n+1)) > 0 end - @testset "Exponential cones" for _ in 1:30 + + @testset "Exponential and power cones" for _ in 1:30 (x, y, z) = rand(3) y += 1 # ensure y > 0 if y * exp(x/y) <= z @@ -82,6 +83,23 @@ import LinearAlgebra elseif u < 0 @test MOI.distance_to_set([u, v, w], MOI.DualExponentialCone()) ≈ u*exp(v/u) - ℯ * w end + + end + (x, y) = randn(2) + if x < 0 || y < 0 + for e in (10 * rand(10) .- 5) # e in [-5, 5] + @test MOI.distance_to_set([x, y, 0.0], MOI.PowerCone(e)) > 0 + end + else + for e in (10 * rand(10) .- 5) # e in [-5, 5] + r = x^e * y^(1-e) + for z in -r:-0.5:r + @test MOI.distance_to_set([x, y, z], MOI.PowerCone(e)) ≈ 0 atol=eps(Float64) + end + @test MOI.distance_to_set([x, y, 3r], MOI.PowerCone(e)) ≈ MOI.distance_to_set([x, y, -3r], MOI.PowerCone(e)) > 0 + end end end + + end From b6f979a7f13fee7d17434e56eb0e9c1f587348fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Fri, 6 Mar 2020 20:12:18 +0100 Subject: [PATCH 22/28] dual power cone --- src/sets.jl | 14 +++++----- test/set_distances.jl | 65 ++++++++++++++++++++++++------------------- 2 files changed, 44 insertions(+), 35 deletions(-) diff --git a/src/sets.jl b/src/sets.jl index 1aa6b5709b..98bfea78ab 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -404,7 +404,7 @@ function distance_to_set(vs::AbstractVector{<:Real}, s::DualExponentialCone) w = vs[3] result = -u*exp(v/u) - ℯ * w return LinearAlgebra.norm2( - (max(-u, zero(result)), max(result, zero(result))) + (max(u, zero(result)), max(result, zero(result))) ) end @@ -429,7 +429,7 @@ function distance_to_set(v::AbstractVector{<:Real}, s::PowerCone) # early return to avoid complex exponent results if x < 0 || y < 0 return LinearAlgebra.norm2( - (max(-x, zero(result)), max(-y, zero(result))) + (max(-x, zero(x)), max(-y, zero(x))) ) end result = abs(z) - x^e * y^(1-e) @@ -450,11 +450,11 @@ 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] +function distance_to_set(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 diff --git a/test/set_distances.jl b/test/set_distances.jl index c1f41dd1b9..eaed44639f 100644 --- a/test/set_distances.jl +++ b/test/set_distances.jl @@ -67,39 +67,48 @@ import LinearAlgebra @test MOI.distance_to_set(vcat(t / 2, -x), MOI.GeometricMeanCone(n+1)) > 0 end - @testset "Exponential and power cones" for _ in 1:30 - (x, y, z) = rand(3) - y += 1 # ensure y > 0 - if y * exp(x/y) <= z - @test MOI.distance_to_set([x, y, z], MOI.ExponentialCone()) ≈ 0 atol=eps(Float64) - @test MOI.distance_to_set([x, -1, z], MOI.ExponentialCone()) ≈ 1 atol=eps(Float64) - else - @test MOI.distance_to_set([x, y, z], MOI.ExponentialCone()) ≈ y * exp(x/y) - z - end - (u, v, w) = randn(3) - if u != 0.0 # just in case not to blow up - if u*exp(v/u) < ℯ * w && u < 0 - @test MOI.distance_to_set([u, v, w], MOI.DualExponentialCone()) ≈ 0 atol=eps(Float64) - elseif u < 0 - @test MOI.distance_to_set([u, v, w], MOI.DualExponentialCone()) ≈ u*exp(v/u) - ℯ * w + @testset "Exponential and power cones" begin + for _ in 1:30 + (x, y, z) = rand(3) + y += 1 # ensure y > 0 + if y * exp(x/y) <= z + @test MOI.distance_to_set([x, y, z], MOI.ExponentialCone()) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set([x, -1, z], MOI.ExponentialCone()) ≈ 1 atol=eps(Float64) + else + @test MOI.distance_to_set([x, y, z], MOI.ExponentialCone()) ≈ y * exp(x/y) - z end - - end - (x, y) = randn(2) - if x < 0 || y < 0 - for e in (10 * rand(10) .- 5) # e in [-5, 5] - @test MOI.distance_to_set([x, y, 0.0], MOI.PowerCone(e)) > 0 + (u, v, w) = randn(3) + if u != 0.0 # just in case not to blow up + if -u*exp(v/u) < ℯ * w && u < 0 + @test MOI.distance_to_set([u, v, w], MOI.DualExponentialCone()) ≈ 0 atol=eps(Float64) + elseif u < 0 + @test MOI.distance_to_set([u, v, w], MOI.DualExponentialCone()) ≈ -u*exp(v/u) - ℯ * w + end + end - else - for e in (10 * rand(10) .- 5) # e in [-5, 5] - r = x^e * y^(1-e) - for z in -r:-0.5:r - @test MOI.distance_to_set([x, y, z], MOI.PowerCone(e)) ≈ 0 atol=eps(Float64) + (x, y) = randn(2) + if x < 0 || y < 0 + for e in (10 * rand(10) .- 5) # e in [-5, 5] + @test MOI.distance_to_set([x, y, 0.0], MOI.PowerCone(e)) > 0 + end + else + for e in (10 * rand(10) .- 5) # e in [-5, 5] + r = x^e * y^(1-e) + for z in -r:-0.5:r + @test MOI.distance_to_set([x, y, z], MOI.PowerCone(e)) ≈ 0 atol=eps(Float64) + end + @test MOI.distance_to_set([x, y, 3r], MOI.PowerCone(e)) ≈ MOI.distance_to_set([x, y, -3r], MOI.PowerCone(e)) > 0 end - @test MOI.distance_to_set([x, y, 3r], MOI.PowerCone(e)) ≈ MOI.distance_to_set([x, y, -3r], MOI.PowerCone(e)) > 0 + end + + (u, v, w) = 10 * rand(3) + e = rand() + if 0 < e < 1 # avoid exponents of negatives + @test MOI.distance_to_set([u, v, 0.0], MOI.DualPowerCone(e)) ≈ 0 atol = 10eps(Float64) + @test MOI.distance_to_set([u, v, u^e * v^(1-e) / (e^e * (1-e)^(1-e))], MOI.DualPowerCone(e)) ≈ 0 atol=20eps(Float64) + @test MOI.distance_to_set([u, v, 1 + u^e * v^(1-e) / (e^e * (1-e)^(1-e))], MOI.DualPowerCone(e)) ≈ 1 end end end - end From b19f884b1e072d994487a5b0d3be745238dbc2c4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Mon, 9 Mar 2020 10:26:12 +0100 Subject: [PATCH 23/28] increased tol --- test/set_distances.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/set_distances.jl b/test/set_distances.jl index eaed44639f..e6103dbba6 100644 --- a/test/set_distances.jl +++ b/test/set_distances.jl @@ -105,7 +105,7 @@ import LinearAlgebra e = rand() if 0 < e < 1 # avoid exponents of negatives @test MOI.distance_to_set([u, v, 0.0], MOI.DualPowerCone(e)) ≈ 0 atol = 10eps(Float64) - @test MOI.distance_to_set([u, v, u^e * v^(1-e) / (e^e * (1-e)^(1-e))], MOI.DualPowerCone(e)) ≈ 0 atol=20eps(Float64) + @test MOI.distance_to_set([u, v, u^e * v^(1-e) / (e^e * (1-e)^(1-e))], MOI.DualPowerCone(e)) ≈ 0 atol=100eps(Float64) @test MOI.distance_to_set([u, v, 1 + u^e * v^(1-e) / (e^e * (1-e)^(1-e))], MOI.DualPowerCone(e)) ≈ 1 end end From 54ff584afd7a5aa378609a0d95ef577e2ccec143 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Sat, 16 May 2020 16:15:49 +0100 Subject: [PATCH 24/28] default distance --- src/sets.jl | 71 +++++++++++++++++++++------------- test/set_distances.jl | 88 +++++++++++++++++++++++++------------------ 2 files changed, 95 insertions(+), 64 deletions(-) diff --git a/src/sets.jl b/src/sets.jl index 86c86a66ea..defbb7c8b7 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -82,14 +82,30 @@ function dual_set_type end dual_set_type(S::Type{<:AbstractSet}) = error("Dual type of $S is not implemented.") + +""" + DistanceFunction + +Distance function used to evaluate the distance from a point to a set. +""" +abstract type DistanceFunction end + """ - distance_to_set(v, s) + DefaultDistance + +Default distance function, uses the Euclidean distance. +""" +struct DefaultDistance <: DistanceFunction end + +""" + distance_to_set(distance_definition, 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. +Each set `S` implements `distance_to_set(d::DistanceFunction, v::T, s::S)` +with `T` of appropriate type for members of the set. """ function distance_to_set end @@ -130,7 +146,7 @@ end dual_set(s::Reals) = Zeros(dimension(s)) dual_set_type(::Type{Reals}) = Zeros -function distance_to_set(v::AbstractVector{T}, s::Reals) where {T <: Real} +function distance_to_set(::DefaultDistance, v::AbstractVector{T}, s::Reals) where {T <: Real} _check_dimension(v, s) return zero(T) end @@ -147,7 +163,7 @@ end dual_set(s::Zeros) = Reals(dimension(s)) dual_set_type(::Type{Zeros}) = Reals -function distance_to_set(v::AbstractVector{<:Real}, s::Zeros) +function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::Zeros) _check_dimension(v, s) return LinearAlgebra.norm2(v) end @@ -164,7 +180,7 @@ end dual_set(s::Nonnegatives) = copy(s) dual_set_type(::Type{Nonnegatives}) = Nonnegatives -function distance_to_set(v::AbstractVector{<:Real}, s::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 @@ -181,7 +197,7 @@ end dual_set(s::Nonpositives) = copy(s) dual_set_type(::Type{Nonpositives}) = Nonpositives -function distance_to_set(v::AbstractVector{<:Real}, s::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 @@ -217,9 +233,9 @@ 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) +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) @@ -248,7 +264,7 @@ 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)) +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}) @@ -271,7 +287,7 @@ end dual_set(s::NormInfinityCone) = NormOneCone(dimension(s)) dual_set_type(::Type{NormInfinityCone}) = NormOneCone -function distance_to_set(v::AbstractVector{<:Real}, s::NormInfinityCone) +function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::NormInfinityCone) _check_dimension(v, s) t = v[1] xs = v[2:end] @@ -291,7 +307,7 @@ end dual_set(s::NormOneCone) = NormInfinityCone(dimension(s)) dual_set_type(::Type{NormOneCone}) = NormInfinityCone -function distance_to_set(v::AbstractVector{<:Real}, s::NormOneCone) +function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::NormOneCone) _check_dimension(v, s) t = v[1] xs = v[2:end] @@ -311,7 +327,7 @@ end dual_set(s::SecondOrderCone) = copy(s) dual_set_type(::Type{SecondOrderCone}) = SecondOrderCone -function distance_to_set(v::AbstractVector{<:Real}, s::SecondOrderCone) +function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::SecondOrderCone) _check_dimension(v, s) t = v[1] xs = v[2:end] @@ -331,7 +347,7 @@ end dual_set(s::RotatedSecondOrderCone) = copy(s) dual_set_type(::Type{RotatedSecondOrderCone}) = RotatedSecondOrderCone -function distance_to_set(v::AbstractVector{<:Real}, s::RotatedSecondOrderCone) +function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::RotatedSecondOrderCone) _check_dimension(v, s) t = v[1] u = v[2] @@ -355,7 +371,7 @@ struct GeometricMeanCone <: AbstractVectorSet dimension::Int end -function distance_to_set(v::AbstractVector{<:Real}, s::GeometricMeanCone) +function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::GeometricMeanCone) _check_dimension(v, s) t = v[1] xs = v[2:end] @@ -381,7 +397,7 @@ struct ExponentialCone <: AbstractVectorSet end dual_set(s::ExponentialCone) = DualExponentialCone() dual_set_type(::Type{ExponentialCone}) = DualExponentialCone -function distance_to_set(v::AbstractVector{<:Real}, s::ExponentialCone) +function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::ExponentialCone) _check_dimension(v, s) x = v[1] y = v[2] @@ -402,7 +418,7 @@ struct DualExponentialCone <: AbstractVectorSet end dual_set(s::DualExponentialCone) = ExponentialCone() dual_set_type(::Type{DualExponentialCone}) = ExponentialCone -function distance_to_set(vs::AbstractVector{<:Real}, s::DualExponentialCone) +function distance_to_set(::DefaultDistance, vs::AbstractVector{<:Real}, s::DualExponentialCone) _check_dimension(vs, s) u = vs[1] v = vs[2] @@ -425,7 +441,7 @@ 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) +function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::PowerCone) _check_dimension(v, s) x = v[1] y = v[2] @@ -455,7 +471,7 @@ 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(vs::AbstractVector{<:Real}, s::DualPowerCone) +function distance_to_set(::DefaultDistance, vs::AbstractVector{<:Real}, s::DualPowerCone) _check_dimension(vs, s) u = vs[1] v = vs[2] @@ -489,7 +505,7 @@ struct RelativeEntropyCone <: AbstractVectorSet dimension::Int end -function distance_to_set(v::AbstractVector{<:Real}, set::RelativeEntropyCone) +function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, set::RelativeEntropyCone) _check_dimension(v, s) n = (dimension(set)-1) ÷ 2 u = v[1] @@ -517,7 +533,7 @@ 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) +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)) @@ -540,7 +556,7 @@ 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) +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)) @@ -832,11 +848,11 @@ The set ``\\{ 0, 1 \\}``. """ struct ZeroOne <: AbstractScalarSet end -function distance_to_set(v::T, ::ZeroOne) where {T <: Real} +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(v::Real, ::Integer) +function distance_to_set(::DefaultDistance, v::Real, ::Integer) return min(abs(v - floor(v)), abs(v - ceil(v))) end @@ -878,7 +894,7 @@ struct SOS1{T <: Real} <: AbstractVectorSet 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} +function distance_to_set(::DefaultDistance, v::AbstractVector{T}, ::SOS1) where {T <: Real} _check_dimension(v, s) d = distance_to_set(v, Zeros(length(v))) m = maximum(d) @@ -960,7 +976,7 @@ end dimension(::IndicatorSet) = 2 # takes in input [z, f(x)] -function distance_to_set(v::AbstractVector{T}, s::IndicatorSet{A}) where {A, T <: Real} +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 @@ -968,7 +984,7 @@ function distance_to_set(v::AbstractVector{T}, s::IndicatorSet{A}) where {A, T < return zeros(T, 2) end return LinearAlgebra.norm2( - (distance_to_set(z, ZeroOne()), distance_to_set(v[2], s.set)) + (distance_to_set(d, z, ZeroOne()), distance_to_set(v[2], s.set)) ) end @@ -1069,6 +1085,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 diff --git a/test/set_distances.jl b/test/set_distances.jl index e6103dbba6..01841be127 100644 --- a/test/set_distances.jl +++ b/test/set_distances.jl @@ -7,64 +7,64 @@ import LinearAlgebra @testset "Set distances" begin @testset "$n-dimensional orthants" for n in 1:3:15 v = rand(n) - @test MOI.distance_to_set(v, MOI.Reals(n)) ≈ 0 atol=eps(Float64) - @test MOI.distance_to_set(v, MOI.Zeros(n)) > 0 - @test MOI.distance_to_set(v, MOI.Nonnegatives(n)) ≈ 0 atol=eps(Float64) - @test MOI.distance_to_set(-v, MOI.Nonpositives(n)) ≈ 0 atol=eps(Float64) - @test MOI.distance_to_set(v, MOI.Nonpositives(n)) ≈ MOI.distance_to_set(-v, MOI.Nonnegatives(n)) > 0 + @test MOI.distance_to_set(MOI.DefaultDistance(), v, MOI.Reals(n)) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), v, MOI.Zeros(n)) > 0 + @test MOI.distance_to_set(MOI.DefaultDistance(), v, MOI.Nonnegatives(n)) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), -v, MOI.Nonpositives(n)) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), v, MOI.Nonpositives(n)) ≈ MOI.distance_to_set(MOI.DefaultDistance(), -v, MOI.Nonnegatives(n)) > 0 end @testset "Scalar comparisons" begin values = rand(10) for v in values - @test MOI.distance_to_set(v, MOI.EqualTo(v)) ≈ 0 atol=eps(Float64) - @test MOI.distance_to_set(-v, MOI.EqualTo(v)) ≈ MOI.distance_to_set(v, MOI.EqualTo(-v)) ≈ 2v - @test MOI.distance_to_set(v, MOI.LessThan(v)) ≈ MOI.distance_to_set(v, MOI.LessThan(v+1)) ≈ 0 - @test MOI.distance_to_set(v, MOI.LessThan(0)) ≈ MOI.distance_to_set(-v, MOI.GreaterThan(0)) ≈ v - @test MOI.distance_to_set(v, MOI.GreaterThan(v)) ≈ MOI.distance_to_set(v+1, MOI.GreaterThan(v+1)) ≈ 0 - @test MOI.distance_to_set(v, MOI.Interval(v,v)) ≈ MOI.distance_to_set(v, MOI.Interval(-v,v)) ≈ 0 - @test MOI.distance_to_set(v, MOI.Interval(-v, 0.0)) ≈ MOI.distance_to_set(-v, MOI.Interval(0.0, v)) ≈ v + @test MOI.distance_to_set(MOI.DefaultDistance(), v, MOI.EqualTo(v)) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), -v, MOI.EqualTo(v)) ≈ MOI.distance_to_set(MOI.DefaultDistance(), v, MOI.EqualTo(-v)) ≈ 2v + @test MOI.distance_to_set(MOI.DefaultDistance(), v, MOI.LessThan(v)) ≈ MOI.distance_to_set(MOI.DefaultDistance(), v, MOI.LessThan(v+1)) ≈ 0 + @test MOI.distance_to_set(MOI.DefaultDistance(), v, MOI.LessThan(0)) ≈ MOI.distance_to_set(MOI.DefaultDistance(), -v, MOI.GreaterThan(0)) ≈ v + @test MOI.distance_to_set(MOI.DefaultDistance(), v, MOI.GreaterThan(v)) ≈ MOI.distance_to_set(MOI.DefaultDistance(), v+1, MOI.GreaterThan(v+1)) ≈ 0 + @test MOI.distance_to_set(MOI.DefaultDistance(), v, MOI.Interval(v,v)) ≈ MOI.distance_to_set(MOI.DefaultDistance(), v, MOI.Interval(-v,v)) ≈ 0 + @test MOI.distance_to_set(MOI.DefaultDistance(), v, MOI.Interval(-v, 0.0)) ≈ MOI.distance_to_set(MOI.DefaultDistance(), -v, MOI.Interval(0.0, v)) ≈ v end end @testset "$n-dimensional norm cones" for n in 2:5:15 x = rand(n) tsum = sum(x) vsum = vcat(tsum, x) - @test MOI.distance_to_set(vsum, MOI.NormOneCone(n+1)) ≈ MOI.distance_to_set(vsum, MOI.NormInfinityCone(n+1)) ≈ 0 + @test MOI.distance_to_set(MOI.DefaultDistance(), vsum, MOI.NormOneCone(n+1)) ≈ MOI.distance_to_set(MOI.DefaultDistance(), vsum, MOI.NormInfinityCone(n+1)) ≈ 0 tmax = maximum(x) vmax = vcat(tmax, x) - @test MOI.distance_to_set(vmax, MOI.NormOneCone(n+1)) > 0 - @test MOI.distance_to_set(vmax, MOI.NormInfinityCone(n+1)) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), vmax, MOI.NormOneCone(n+1)) > 0 + @test MOI.distance_to_set(MOI.DefaultDistance(), vmax, MOI.NormInfinityCone(n+1)) ≈ 0 atol=eps(Float64) tmin = 0 vmin = vcat(tmin, x) - @test MOI.distance_to_set(vmin, MOI.NormInfinityCone(n+1)) ≈ tmax - @test MOI.distance_to_set(vmin, MOI.NormOneCone(n+1)) ≈ tsum + @test MOI.distance_to_set(MOI.DefaultDistance(), vmin, MOI.NormInfinityCone(n+1)) ≈ tmax + @test MOI.distance_to_set(MOI.DefaultDistance(), vmin, MOI.NormOneCone(n+1)) ≈ tsum tvalid = sqrt(n) # upper bound on the norm2 vok_soc = vcat(tvalid, x) - @test MOI.distance_to_set(vok_soc, MOI.SecondOrderCone(n+1) ) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), vok_soc, MOI.SecondOrderCone(n+1) ) ≈ 0 atol=eps(Float64) vko_soc = vcat(-2, x) - @test MOI.distance_to_set(vko_soc, MOI.SecondOrderCone(n+1) ) ≈ 2 + LinearAlgebra.norm2(x) + @test MOI.distance_to_set(MOI.DefaultDistance(), vko_soc, MOI.SecondOrderCone(n+1) ) ≈ 2 + LinearAlgebra.norm2(x) vko_soc = vcat(-2, x) - @test MOI.distance_to_set(vko_soc, MOI.SecondOrderCone(n+1) ) ≈ 2 + LinearAlgebra.norm2(x) + @test MOI.distance_to_set(MOI.DefaultDistance(), vko_soc, MOI.SecondOrderCone(n+1) ) ≈ 2 + LinearAlgebra.norm2(x) t_ko_rot = u_ko_rot = LinearAlgebra.norm2(x) / 2 vko_roc = vcat(t_ko_rot, u_ko_rot, x) - @test MOI.distance_to_set(vko_roc, MOI.RotatedSecondOrderCone(n+2)) ≈ LinearAlgebra.dot(x,x) / 2 + @test MOI.distance_to_set(MOI.DefaultDistance(), vko_roc, MOI.RotatedSecondOrderCone(n+2)) ≈ LinearAlgebra.dot(x,x) / 2 vok_roc = vcat(t_ko_rot * 2, u_ko_rot, x) - @test MOI.distance_to_set(vok_roc, MOI.RotatedSecondOrderCone(n+2)) ≈ 0 atol=5eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), vok_roc, MOI.RotatedSecondOrderCone(n+2)) ≈ 0 atol=5eps(Float64) end @testset "Geometric Mean cone dimension $n" for n in 2:5:15 x = rand(n) t = 0.5 * prod(x)^(inv(n)) vok = vcat(t, x) - @test MOI.distance_to_set(vok, MOI.GeometricMeanCone(n+1)) ≈ 0 atol=eps(Float64) - @test MOI.distance_to_set(vcat(t / 2, x), MOI.GeometricMeanCone(n+1)) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), vok, MOI.GeometricMeanCone(n+1)) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), vcat(t / 2, x), MOI.GeometricMeanCone(n+1)) ≈ 0 atol=eps(Float64) # negative x always means positive distance - @test MOI.distance_to_set(vcat(t / 2, vcat(x, -1)), MOI.GeometricMeanCone(n+2)) > 0 - @test MOI.distance_to_set(vcat(t / 2, -x), MOI.GeometricMeanCone(n+1)) > 0 + @test MOI.distance_to_set(MOI.DefaultDistance(), vcat(t / 2, vcat(x, -1)), MOI.GeometricMeanCone(n+2)) > 0 + @test MOI.distance_to_set(MOI.DefaultDistance(), vcat(t / 2, -x), MOI.GeometricMeanCone(n+1)) > 0 end @testset "Exponential and power cones" begin @@ -72,43 +72,57 @@ import LinearAlgebra (x, y, z) = rand(3) y += 1 # ensure y > 0 if y * exp(x/y) <= z - @test MOI.distance_to_set([x, y, z], MOI.ExponentialCone()) ≈ 0 atol=eps(Float64) - @test MOI.distance_to_set([x, -1, z], MOI.ExponentialCone()) ≈ 1 atol=eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), [x, y, z], MOI.ExponentialCone()) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), [x, -1, z], MOI.ExponentialCone()) ≈ 1 atol=eps(Float64) else - @test MOI.distance_to_set([x, y, z], MOI.ExponentialCone()) ≈ y * exp(x/y) - z + @test MOI.distance_to_set(MOI.DefaultDistance(), [x, y, z], MOI.ExponentialCone()) ≈ y * exp(x/y) - z end (u, v, w) = randn(3) if u != 0.0 # just in case not to blow up if -u*exp(v/u) < ℯ * w && u < 0 - @test MOI.distance_to_set([u, v, w], MOI.DualExponentialCone()) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), [u, v, w], MOI.DualExponentialCone()) ≈ 0 atol=eps(Float64) elseif u < 0 - @test MOI.distance_to_set([u, v, w], MOI.DualExponentialCone()) ≈ -u*exp(v/u) - ℯ * w + @test MOI.distance_to_set(MOI.DefaultDistance(), [u, v, w], MOI.DualExponentialCone()) ≈ -u*exp(v/u) - ℯ * w end end (x, y) = randn(2) if x < 0 || y < 0 for e in (10 * rand(10) .- 5) # e in [-5, 5] - @test MOI.distance_to_set([x, y, 0.0], MOI.PowerCone(e)) > 0 + @test MOI.distance_to_set(MOI.DefaultDistance(), [x, y, 0.0], MOI.PowerCone(e)) > 0 end else for e in (10 * rand(10) .- 5) # e in [-5, 5] r = x^e * y^(1-e) for z in -r:-0.5:r - @test MOI.distance_to_set([x, y, z], MOI.PowerCone(e)) ≈ 0 atol=eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), [x, y, z], MOI.PowerCone(e)) ≈ 0 atol=eps(Float64) end - @test MOI.distance_to_set([x, y, 3r], MOI.PowerCone(e)) ≈ MOI.distance_to_set([x, y, -3r], MOI.PowerCone(e)) > 0 + @test MOI.distance_to_set(MOI.DefaultDistance(), [x, y, 3r], MOI.PowerCone(e)) ≈ MOI.distance_to_set(MOI.DefaultDistance(), [x, y, -3r], MOI.PowerCone(e)) > 0 end end (u, v, w) = 10 * rand(3) e = rand() if 0 < e < 1 # avoid exponents of negatives - @test MOI.distance_to_set([u, v, 0.0], MOI.DualPowerCone(e)) ≈ 0 atol = 10eps(Float64) - @test MOI.distance_to_set([u, v, u^e * v^(1-e) / (e^e * (1-e)^(1-e))], MOI.DualPowerCone(e)) ≈ 0 atol=100eps(Float64) - @test MOI.distance_to_set([u, v, 1 + u^e * v^(1-e) / (e^e * (1-e)^(1-e))], MOI.DualPowerCone(e)) ≈ 1 + @test MOI.distance_to_set(MOI.DefaultDistance(), [u, v, 0.0], MOI.DualPowerCone(e)) ≈ 0 atol = 10eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), [u, v, u^e * v^(1-e) / (e^e * (1-e)^(1-e))], MOI.DualPowerCone(e)) ≈ 0 atol=100eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), [u, v, 1 + u^e * v^(1-e) / (e^e * (1-e)^(1-e))], MOI.DualPowerCone(e)) ≈ 1 end end end end + +struct DummyDistance <: MOI.DistanceFunction end + +MOI.distance_to_set(::DummyDistance, v, s) = MOI.distance_to_set(MOI.DefaultDistance(), v, s) / 2 + +@testset "Set non-default distance" begin + for n in 1:3 + v = rand(n) + for s in (MOI.Reals(n), MOI.Zeros(n), MOI.Nonnegatives(n), MOI.Nonpositives(n)) + @test MOI.distance_to_set(MOI.DefaultDistance(), v, MOI.Reals(n)) ≈ 2 * MOI.distance_to_set(DummyDistance(), v, MOI.Reals(n)) atol=eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), -v, MOI.Reals(n)) ≈ 2 * MOI.distance_to_set(DummyDistance(), -v, MOI.Reals(n)) atol=eps(Float64) + end + end +end From 68888b5662c4d09033fd059af9fd34fe6c2830af Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Sat, 16 May 2020 17:05:35 +0100 Subject: [PATCH 25/28] abstract and default distance for sets, added complementarity and dummy test distance --- src/sets.jl | 25 +++++++++++++++++++++---- test/set_distances.jl | 18 +++++++++++++++++- 2 files changed, 38 insertions(+), 5 deletions(-) diff --git a/src/sets.jl b/src/sets.jl index defbb7c8b7..5117f2feef 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -84,18 +84,20 @@ dual_set_type(S::Type{<:AbstractSet}) = error("Dual type of $S is not implemente """ - DistanceFunction + 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 DistanceFunction end +abstract type AbstractDistance end """ DefaultDistance Default distance function, uses the Euclidean distance. """ -struct DefaultDistance <: DistanceFunction end +struct DefaultDistance <: AbstractDistance end """ distance_to_set(distance_definition, v, s) @@ -104,7 +106,7 @@ 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(d::DistanceFunction, v::T, s::S)` +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 @@ -1051,6 +1053,21 @@ end dimension(set::Complements) = 2 * set.dimension +# warning, only works for classical Complementarity constraints: +# 0 <= x_i ⟂ F_i(x) >= 0 +function distance_to_set(d::DefaultDistance, v, s::Complements) + _check_dimension(v, s) + non_positives_var = [distance_to_set(d, v[i], GreaterThan(0.0)) for i in 1:s.dimension] + non_positives_func = [distance_to_set(d, v[s.dimension + i], GreaterThan(0.0)) for i in 1:s.dimension] + comp_distance = [ + min(distance_to_set(d, v[i], EqualTo(0.0)), distance_to_set(d, v[i + s.dimension], EqualTo(0.0))) + for i in 1:s.dimension + ] + return LinearAlgebra.norm2( + comp_distance .+ non_positives_func .+ non_positives_var + ) +end + # isbits types, nothing to copy function Base.copy( set::Union{ diff --git a/test/set_distances.jl b/test/set_distances.jl index 01841be127..8f83938ee5 100644 --- a/test/set_distances.jl +++ b/test/set_distances.jl @@ -110,10 +110,26 @@ import LinearAlgebra end end end + + @testset "Complements $n" for n in 2:10 + xs = rand(n) + fs = rand(n) + vals_wrong = vcat(xs, fs) + @test MOI.distance_to_set(MOI.DefaultDistance(), vals_wrong, MOI.Complements(n)) > 0 + xs .*= 0 + vals_trivial_x = vcat(xs, fs) + @test MOI.distance_to_set(MOI.DefaultDistance(), vals_trivial_x, MOI.Complements(n)) ≈ 0 + xs = fs + fs .*= 0 + vals_trivial_f = vcat(xs, fs) + @test MOI.distance_to_set(MOI.DefaultDistance(), vals_trivial_f, MOI.Complements(n)) ≈ 0 + vals_trivial_f[1] = -1 + @test MOI.distance_to_set(MOI.DefaultDistance(), vals_trivial_f, MOI.Complements(n)) > 0 + end end -struct DummyDistance <: MOI.DistanceFunction end +struct DummyDistance <: MOI.AbstractDistance end MOI.distance_to_set(::DummyDistance, v, s) = MOI.distance_to_set(MOI.DefaultDistance(), v, s) / 2 From 21c759f926a299bf791c93586b0345688cc52327 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Sat, 16 May 2020 21:10:44 +0100 Subject: [PATCH 26/28] remove complementarity --- src/sets.jl | 15 --------------- test/set_distances.jl | 19 +------------------ 2 files changed, 1 insertion(+), 33 deletions(-) diff --git a/src/sets.jl b/src/sets.jl index 5117f2feef..36cb7130bf 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -1053,21 +1053,6 @@ end dimension(set::Complements) = 2 * set.dimension -# warning, only works for classical Complementarity constraints: -# 0 <= x_i ⟂ F_i(x) >= 0 -function distance_to_set(d::DefaultDistance, v, s::Complements) - _check_dimension(v, s) - non_positives_var = [distance_to_set(d, v[i], GreaterThan(0.0)) for i in 1:s.dimension] - non_positives_func = [distance_to_set(d, v[s.dimension + i], GreaterThan(0.0)) for i in 1:s.dimension] - comp_distance = [ - min(distance_to_set(d, v[i], EqualTo(0.0)), distance_to_set(d, v[i + s.dimension], EqualTo(0.0))) - for i in 1:s.dimension - ] - return LinearAlgebra.norm2( - comp_distance .+ non_positives_func .+ non_positives_var - ) -end - # isbits types, nothing to copy function Base.copy( set::Union{ diff --git a/test/set_distances.jl b/test/set_distances.jl index 8f83938ee5..30635a48f6 100644 --- a/test/set_distances.jl +++ b/test/set_distances.jl @@ -109,24 +109,7 @@ import LinearAlgebra @test MOI.distance_to_set(MOI.DefaultDistance(), [u, v, 1 + u^e * v^(1-e) / (e^e * (1-e)^(1-e))], MOI.DualPowerCone(e)) ≈ 1 end end - end - - @testset "Complements $n" for n in 2:10 - xs = rand(n) - fs = rand(n) - vals_wrong = vcat(xs, fs) - @test MOI.distance_to_set(MOI.DefaultDistance(), vals_wrong, MOI.Complements(n)) > 0 - xs .*= 0 - vals_trivial_x = vcat(xs, fs) - @test MOI.distance_to_set(MOI.DefaultDistance(), vals_trivial_x, MOI.Complements(n)) ≈ 0 - xs = fs - fs .*= 0 - vals_trivial_f = vcat(xs, fs) - @test MOI.distance_to_set(MOI.DefaultDistance(), vals_trivial_f, MOI.Complements(n)) ≈ 0 - vals_trivial_f[1] = -1 - @test MOI.distance_to_set(MOI.DefaultDistance(), vals_trivial_f, MOI.Complements(n)) > 0 - end - + end end struct DummyDistance <: MOI.AbstractDistance end From 9a2e898d0fe4cf29e8e5503efa9be63ae2ca25f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Tue, 19 May 2020 07:33:45 +0100 Subject: [PATCH 27/28] increase tolerance --- test/set_distances.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/set_distances.jl b/test/set_distances.jl index 30635a48f6..aeee354581 100644 --- a/test/set_distances.jl +++ b/test/set_distances.jl @@ -53,7 +53,7 @@ import LinearAlgebra vko_roc = vcat(t_ko_rot, u_ko_rot, x) @test MOI.distance_to_set(MOI.DefaultDistance(), vko_roc, MOI.RotatedSecondOrderCone(n+2)) ≈ LinearAlgebra.dot(x,x) / 2 vok_roc = vcat(t_ko_rot * 2, u_ko_rot, x) - @test MOI.distance_to_set(MOI.DefaultDistance(), vok_roc, MOI.RotatedSecondOrderCone(n+2)) ≈ 0 atol=5eps(Float64) + @test MOI.distance_to_set(MOI.DefaultDistance(), vok_roc, MOI.RotatedSecondOrderCone(n+2)) ≈ 0 atol=10eps(Float64) end @testset "Geometric Mean cone dimension $n" for n in 2:5:15 From 13725c30e340d3d2aa21390e59cea5922f62929d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mathieu=20Besan=C3=A7on?= Date: Tue, 19 May 2020 16:17:08 +0100 Subject: [PATCH 28/28] update doc --- src/sets.jl | 22 ++++++---------------- 1 file changed, 6 insertions(+), 16 deletions(-) diff --git a/src/sets.jl b/src/sets.jl index 36cb7130bf..aa5089a3ff 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -2,7 +2,6 @@ # Note: When adding a new set, also add it to Utilities.Model. import LinearAlgebra -using LinearAlgebra: dot """ AbstractSet @@ -103,7 +102,6 @@ struct DefaultDistance <: AbstractDistance end distance_to_set(distance_definition, 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 at least `distance_to_set(d::DefaultDistance, v::T, s::S)` @@ -111,6 +109,8 @@ 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 @@ -334,7 +334,7 @@ function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::Second t = v[1] xs = v[2:end] result = LinearAlgebra.norm2(xs) - t - return max(result, zero(result)) # avoids sqrt + return max(result, zero(result)) end """ @@ -355,7 +355,7 @@ function distance_to_set(::DefaultDistance, v::AbstractVector{<:Real}, s::Rotate u = v[2] xs = v[3:end] return LinearAlgebra.norm2( - (max(-t, zero(t)), max(-u, zero(u)), max(dot(xs,xs) - 2 * t * u)) + (max(-t, zero(t)), max(-u, zero(u)), max(LinearAlgebra.dot(xs,xs) - 2 * t * u)) ) end @@ -898,18 +898,8 @@ 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) - 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 + _, i = findmax(abs.(v)) + return LinearAlgebra.norm2([v[j] for j = eachindex(v) if j != i]) end """