Merge pull request #9 from JuliaAlgebra/bl/symmetric_setindex!

Add symmetric_setindex!

docs/src/atoms.md

@@ -5,6 +5,8 @@
 ```@docs
 MomentMatrix
 moment_matrix
+SymMatrix
+symmetric_setindex!
 ```
 
 ## Atomic measure

src/MultivariateMoments.jl

@@ -16,6 +16,7 @@ abstract type AbstractMeasure{T} <: AbstractMeasureLike{T} end
 include("moment.jl")
 include("measure.jl")
 include("expectation.jl")
+include("symmatrix.jl")
 include("moment_matrix.jl")
 include("atomic.jl")
 include("extract.jl")

src/moment_matrix.jl

@@ -1,40 +1,7 @@
-export SymMatrix, MomentMatrix, getmat, moment_matrix
+export SymMatrix, MomentMatrix, getmat, moment_matrix, symmetric_setindex!
 
 using SemialgebraicSets
 
-struct SymMatrix{T} <: AbstractMatrix{T}
-    Q::Vector{T}
-    n::Int
-end
-
-# i <= j
-#function trimap(i, j, n) # It was used when Q was the lower triangular part
-#    div(n*(n+1), 2) - div((n-i+1)*(n-i+2), 2) + j-i+1
-#end
-
-# j <= i
-function trimap(i, j)
-    div((i-1)*i, 2) + j
-end
-
-function trimat(::Type{T}, f, n, σ) where {T}
-    Q = Vector{T}(undef, trimap(n, n))
-    for i in 1:n
-        for j in 1:i
-            Q[trimap(i, j)] = f(σ[i], σ[j])
-        end
-    end
-    SymMatrix{T}(Q, n)
-end
-
-Base.size(Q::SymMatrix) = (Q.n, Q.n)
-
-function Base.getindex(Q::SymMatrix, i::Integer, j::Integer)
-    Q.Q[trimap(max(i, j), min(i, j))]
-end
-Base.getindex(Q::SymMatrix, I::Tuple) = Q[I...]
-Base.getindex(Q::SymMatrix, I::CartesianIndex) = Q[I.I]
-
 """
     mutable struct MomentMatrix{T, MT <: AbstractMonomial, MVT <: AbstractVector{MT}} <: AbstractMeasureLike{T}
         Q::SymMatrix{T}

src/symmatrix.jl

@@ -0,0 +1,54 @@
+"""
+    struct SymMatrix{T} <: AbstractMatrix{T}
+        Q::Vector{T}
+        n::Int
+    end
+
+Symmetric ``n \\times n`` matrix storing the vectorized upper triangular
+part of the matrix in the `Q` vector (this corresponds to the vectorized
+form of `MathOptInterface.PositiveSemidefiniteConeTriangle`).
+It implement the `AbstractMatrix` interface except for `setindex!` as it might
+break its symmetry. The [`symmetric_setindex!`](@ref) function should be used
+instead.
+"""
+struct SymMatrix{T} <: AbstractMatrix{T}
+    Q::Vector{T}
+    n::Int
+end
+
+# i <= j
+#function trimap(i, j, n) # It was used when Q was the lower triangular part
+#    div(n*(n+1), 2) - div((n-i+1)*(n-i+2), 2) + j-i+1
+#end
+
+# j <= i
+function trimap(i, j)
+    div((i-1)*i, 2) + j
+end
+
+function trimat(::Type{T}, f, n, σ) where {T}
+    Q = Vector{T}(undef, trimap(n, n))
+    for i in 1:n
+        for j in 1:i
+            Q[trimap(i, j)] = f(σ[i], σ[j])
+        end
+    end
+    SymMatrix{T}(Q, n)
+end
+
+Base.size(Q::SymMatrix) = (Q.n, Q.n)
+
+"""
+    symmetric_setindex!(Q::SymMatrix, value, i::Integer, j::Integer)
+
+Set `Q[i, j]` and `Q[j, i]` to the value `value`.
+"""
+function symmetric_setindex!(Q::SymMatrix, value, i::Integer, j::Integer)
+    Q.Q[trimap(max(i, j), min(i, j))] = value
+end
+
+function Base.getindex(Q::SymMatrix, i::Integer, j::Integer)
+    return Q.Q[trimap(max(i, j), min(i, j))]
+end
+Base.getindex(Q::SymMatrix, I::Tuple) = Q[I...]
+Base.getindex(Q::SymMatrix, I::CartesianIndex) = Q[I.I]

test/moment_matrix.jl

@@ -1,3 +1,13 @@
+@testset "SymMatrix" begin
+    Q = SymMatrix([1, 2, 3], 2)
+    symmetric_setindex!(Q, 4, 1, 1)
+    @test Q.Q == [4, 2, 3]
+    symmetric_setindex!(Q, 5, 1, 2)
+    @test Q.Q == [4, 5, 3]
+    symmetric_setindex!(Q, 6, 2, 2)
+    @test Q.Q == [4, 5, 6]
+end
+
 @testset "MomentMatrix" begin
     Mod.@polyvar x y
     @test_throws ArgumentError moment_matrix(measure([1], [x]), [y])