From 39e31f265d2aa7d7c90c9e84c6c1e083d4cbfe03 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 6 Jul 2023 18:49:19 +0200 Subject: [PATCH] Fixes --- src/Utilities/lazy_iterators.jl | 4 ++ src/Utilities/results.jl | 3 +- src/Utilities/sets.jl | 65 +++++++++++++++++++++++++++++++-- test/Utilities/sets.jl | 11 +++++- 4 files changed, 77 insertions(+), 6 deletions(-) diff --git a/src/Utilities/lazy_iterators.jl b/src/Utilities/lazy_iterators.jl index 2d434e1788..558f2d3260 100644 --- a/src/Utilities/lazy_iterators.jl +++ b/src/Utilities/lazy_iterators.jl @@ -53,3 +53,7 @@ Base.eltype(::LazyMap{T}) where {T} = T function Iterators.reverse(it::LazyMap{T}) where {T} return LazyMap{T}(it.f, Iterators.reverse(it.data)) end + +function Base.getindex(it::LazyMap, i) + return it.f(getindex(it.data, i)) +end diff --git a/src/Utilities/results.jl b/src/Utilities/results.jl index b56820d9b1..c6c8aa7bd2 100644 --- a/src/Utilities/results.jl +++ b/src/Utilities/results.jl @@ -547,7 +547,8 @@ function MOI.Utilities.set_dot( set::MOI.HermitianPositiveSemidefiniteConeTriangle, ) sym = MOI.PositiveSemidefiniteConeTriangle(set.side_dimension) - result = set_dot(x, y, sym) + I = Base.OneTo(MOI.dimension(sym)) + result = set_dot(view(x, I), view(y, I), sym) for k in (MOI.dimension(sym)+1):MOI.dimension(set) result = MA.add_mul!!(result, 2, x[k], y[k]) end diff --git a/src/Utilities/sets.jl b/src/Utilities/sets.jl index 8b2e702b37..980502c6f0 100644 --- a/src/Utilities/sets.jl +++ b/src/Utilities/sets.jl @@ -158,6 +158,20 @@ function trimap(row::Integer, column::Integer) return div((row - 1) * row, 2) + column end +struct ZeroVector{T} <: AbstractVector{T} + n::Int +end + +Base.eltype(::Type{ZeroVector{T}}) where {T} = T + +Base.length(v::ZeroVector) = v.n + +Base.size(v::ZeroVector) = (v.n,) + +function Base.getindex(v::ZeroVector{T}, i::Integer) where {T} + return zero(T) +end + struct CanonicalVector{T} <: AbstractVector{T} index::Int n::Int @@ -173,6 +187,14 @@ function Base.getindex(v::CanonicalVector{T}, i::Integer) where {T} return convert(T, i == v.index) end +function Base.view(v::CanonicalVector{T}, I::AbstractUnitRange) where {T} + if v.index in I + return CanonicalVector{T}(v.index - first(I) + 1, length(I)) + else + return ZeroVector{T}(length(I)) + end +end + # This is much faster than the default implementation that goes # through all entries even if only one is nonzero. function LinearAlgebra.dot( @@ -220,6 +242,9 @@ struct SetDotScalingVector{T,S<:MOI.AbstractSet} <: AbstractVector{T} set::S len::Int end +function SetDotScalingVector{T}(s::MOI.AbstractSet, len) where {T} + return SetDotScalingVector{T,typeof(s)}(s, len) +end function Base.getindex(s::SetDotScalingVector{T}, i::Base.Integer) where {T} return sqrt(_set_dot(i, s.set, T)) @@ -228,15 +253,49 @@ end Base.size(x::SetDotScalingVector) = (x.len,) function symmetric_matrix_scaling_vector(::Type{T}, n) where {T} - d = MOI.side_dimension_for_vectorized_dimension(n) - set = MOI.ScaledPositiveSemidefiniteConeTriangle(d) - return SetDotScalingVector(set, n) + d = side_dimension_for_vectorized_dimension(n) + set = MOI.PositiveSemidefiniteConeTriangle(d) + return SetDotScalingVector{T}(set, n) end function symmetric_matrix_inverse_scaling_vector(::Type{T}, n) where {T} return LazyMap{T}(inv, symmetric_matrix_scaling_vector(T, n)) end +""" + struct SymmetricMatrixScalingVector{T} <: AbstractVector{T} + no_scaling::T + scaling::T + len::Int + end + +Vector of scaling for the entries of the vectorized form of +a symmetric matrix. The values `no_scaling` and `scaling` +are stored in the `struct` to avoid creating a new one for each entry. + +!!! warning + This type is deprecated, use `SetDotScalingVector` instead. +""" +struct SymmetricMatrixScalingVector{T} <: AbstractVector{T} + scaling::T + no_scaling::T + len::Int +end + +function SymmetricMatrixScalingVector{T}(scaling::T, len::Int) where {T} + return SymmetricMatrixScalingVector{T}(scaling, one(T), len) +end + +function Base.getindex(s::SymmetricMatrixScalingVector, i::Base.Integer) + if is_diagonal_vectorized_index(i) + return s.no_scaling + else + return s.scaling + end +end + +Base.size(x::SymmetricMatrixScalingVector) = (x.len,) + similar_type(::Type{<:MOI.LessThan}, ::Type{T}) where {T} = MOI.LessThan{T} function similar_type(::Type{<:MOI.GreaterThan}, ::Type{T}) where {T} diff --git a/test/Utilities/sets.jl b/test/Utilities/sets.jl index e2e85218dd..cf8be7c35e 100644 --- a/test/Utilities/sets.jl +++ b/test/Utilities/sets.jl @@ -197,23 +197,29 @@ function test_trimap() return end -function test_symmetric_matrix_scaling() - n = 10 +function test_set_dot_scaling(n = 10) N = div(n * (n + 1), 2) + M = N + div(n * (n - 1), 2) v = MOI.Utilities.SymmetricMatrixScalingVector{Float64}(1.5, 0.5, N) s = MOI.Utilities.symmetric_matrix_scaling_vector(Float64, N) s32 = MOI.Utilities.symmetric_matrix_scaling_vector(Float32, N) is = MOI.Utilities.symmetric_matrix_inverse_scaling_vector(Float64, N) is32 = MOI.Utilities.symmetric_matrix_inverse_scaling_vector(Float32, N) + hpsd = MOI.HermitianPositiveSemidefiniteConeTriangle(n) + hermitian = MOI.Utilities.SetDotScalingVector{Float64}(hpsd, N) k = 0 + imag_k = 0 for j in 1:n for i in 1:(j-1) k += 1 + imag_k += 1 @test v[k] == 1.5 @test s[k] == √2 @test s32[k] == √Float32(2) @test is[k] == inv(√2) @test is32[k] == inv(√Float32(2)) + @test hermitian[k] == √2 + @test hermitian[N+imag_k] == √2 end k += 1 @test v[k] == 0.5 @@ -221,6 +227,7 @@ function test_symmetric_matrix_scaling() @test isone(s32[k]) @test isone(is[k]) @test isone(is32[k]) + @test isone(hermitian[k]) end return end