From 7f212f6ea5c7de993c20422173f6fb6b027241c4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 20 Jan 2020 15:11:15 +0100 Subject: [PATCH 1/2] Add support for arbitrary polynomial basis --- .travis.yml | 4 ++- Project.toml | 3 +- src/MultivariateMoments.jl | 3 ++ src/extract.jl | 15 +++++----- src/moment_matrix.jl | 57 +++++++++++++++++++------------------- test/moment_matrix.jl | 6 ++-- test/runtests.jl | 4 ++- 7 files changed, 51 insertions(+), 41 deletions(-) diff --git a/.travis.yml b/.travis.yml index ec9797d..c81d48e 100644 --- a/.travis.yml +++ b/.travis.yml @@ -5,11 +5,13 @@ os: - osx julia: - 1.0 - - 1.1 + - 1.3 notifications: email: false git: depth: 99999999 +before_script: + - julia -e 'using Pkg; Pkg.add(PackageSpec(url="https://github.com/JuliaAlgebra/MultivariateBases.jl", rev="master"))' after_success: - julia -e 'using Pkg; Pkg.add("Coverage"); cd(Pkg.dir("MultivariateMoments")); using Coverage; Coveralls.submit(process_folder()); Codecov.submit(process_folder())' jobs: diff --git a/Project.toml b/Project.toml index d675734..07bfa7f 100644 --- a/Project.toml +++ b/Project.toml @@ -5,14 +5,15 @@ version = "0.2.4" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378" MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" RowEchelon = "af85af4c-bcd5-5d23-b03a-a909639aa875" SemialgebraicSets = "8e049039-38e8-557d-ae3a-bc521ccf6204" [compat] MultivariatePolynomials = "0.3" -SemialgebraicSets = "0.2" RowEchelon = "0.1" +SemialgebraicSets = "0.2" julia = "1" [extras] diff --git a/src/MultivariateMoments.jl b/src/MultivariateMoments.jl index 0a1f1fc..b81158c 100644 --- a/src/MultivariateMoments.jl +++ b/src/MultivariateMoments.jl @@ -5,6 +5,9 @@ using LinearAlgebra using MultivariatePolynomials const MP = MultivariatePolynomials +import MultivariateBases +const MB = MultivariateBases + export AbstractMeasureLike, AbstractMomentLike abstract type AbstractMeasureLike{T} end abstract type AbstractMomentLike{T} <: AbstractMeasureLike{T} end diff --git a/src/extract.jl b/src/extract.jl index 9e96f01..aed46e8 100644 --- a/src/extract.jl +++ b/src/extract.jl @@ -45,13 +45,13 @@ using SemialgebraicSets # r, vals #end -function build_system(U::AbstractMatrix, mv::AbstractVector, ztol) +function build_system(U::AbstractMatrix, basis::MB.MonomialBasis, ztol) # System is # y = [U 0] * y # where y = x[end:-1:1] # which is # y = U * β - m = length(mv) + m = length(basis) r = size(U, 1) pivots = [findfirst(j -> U[i, j] != 0, 1:m) for i in 1:r] if any(pivots .== 0) @@ -60,7 +60,8 @@ function build_system(U::AbstractMatrix, mv::AbstractVector, ztol) r = length(pivots) U = U[keep, :] end - β = monovec(mv[m + 1 .- pivots]) # monovec makes sure it stays sorted, TypedPolynomials wouldn't let it sorted + monos = basis.monomials + β = monovec(monos[m + 1 .- pivots]) # monovec makes sure it stays sorted, TypedPolynomials wouldn't let it sorted function equation(i) if iszero(r) # sum throws ArgumentError: reducing over an empty collection is not allowed, if r is zero z = zero(eltype(β)) * zero(eltype(U)) @@ -68,11 +69,11 @@ function build_system(U::AbstractMatrix, mv::AbstractVector, ztol) else s = sum(j -> β[r+1-j] * U[j, i], 1:r) end - s - mv[m+1-i] + s - monos[m+1-i] end - system = filter(p -> maxdegree(p) > 0, map(equation, 1:length(mv))) + system = filter(p -> maxdegree(p) > 0, map(equation, 1:length(monos))) # Type instability here :( - if mindegree(mv) == maxdegree(mv) # Homogeneous + if mindegree(monos) == maxdegree(monos) # Homogeneous projectivealgebraicset(system, Buchberger(ztol)) else algebraicset(system, Buchberger(ztol)) @@ -158,7 +159,7 @@ function computesupport!(μ::MomentMatrix, ranktol::Real, dec::LowRankChol=SVDCh # so we take √||M|| = √nM rref!(W, √(nM) * cM / sqrt(m)) #r, vals = solve_system(U', μ.x) - μ.support = build_system(W, μ.x, √cM) # TODO determine what is better between ranktol and sqrt(ranktol) here + μ.support = build_system(W, μ.basis, √cM) # TODO determine what is better between ranktol and sqrt(ranktol) here end """ diff --git a/src/moment_matrix.jl b/src/moment_matrix.jl index 64009d4..ce68b7c 100644 --- a/src/moment_matrix.jl +++ b/src/moment_matrix.jl @@ -3,38 +3,38 @@ export getmat, moment_matrix, symmetric_setindex! using SemialgebraicSets -abstract type AbstractMomentMatrix{T, MT <: AbstractMonomial, MVT <: AbstractVector{MT}} <: AbstractMeasureLike{T} end +abstract type AbstractMomentMatrix{T, B<:MB.AbstractPolynomialBasis} <: AbstractMeasureLike{T} end """ - mutable struct MomentMatrix{T, MT <: AbstractMonomial, MVT <: AbstractVector{MT}} <: AbstractMeasureLike{T} + mutable struct MomentMatrix{T, B<:MultivariateBases.AbstractPolynomialBasis} <: AbstractMeasureLike{T} Q::SymMatrix{T} - x::MVT + basis::B support::Union{Nothing, AlgebraicSet} end Measure ``\\nu`` represented by the moments of the monomial matrix ``x x^\\top`` in the symmetric matrix `Q`. The set of points that are zeros of all the polynomials ``p`` such that ``\\mathbb{E}_{\\nu}[p] = 0`` is stored in the field `support` when it is computed. """ -mutable struct MomentMatrix{T, MT <: AbstractMonomial, MVT <: AbstractVector{MT}} <: AbstractMomentMatrix{T, MT, MVT} +mutable struct MomentMatrix{T, B<:MB.AbstractPolynomialBasis} <: AbstractMomentMatrix{T, B} Q::SymMatrix{T} - x::MVT + basis::B support::Union{Nothing, AlgebraicSet} end -MomentMatrix{T, MT, MVT}(Q::SymMatrix{T}, x::MVT) where {T, MT, MVT} = MomentMatrix{T, MT, MVT}(Q, x, nothing) +MomentMatrix{T, B}(Q::SymMatrix{T}, basis::MB.AbstractPolynomialBasis) where {T, B} = MomentMatrix{T, B}(Q, basis, nothing) +function MomentMatrix(Q::SymMatrix{T}, basis::MB.AbstractPolynomialBasis) where T + return MomentMatrix{T, typeof(basis)}(Q, basis) +end -MP.variables(μ::MomentMatrix) = variables(μ.x) -MP.nvariables(μ::MomentMatrix) = nvariables(μ.x) +MP.variables(μ::MomentMatrix) = variables(μ.basis) +MP.nvariables(μ::MomentMatrix) = nvariables(μ.basis) -function MomentMatrix{T}(f::Function, x::AbstractVector{MT}, σ) where {T, MT<:AbstractMonomial} - MomentMatrix{T, MT, monovectype(x)}(trimat(T, f, length(x), σ), x) +function MomentMatrix{T}(f::Function, basis::MB.AbstractPolynomialBasis, σ=1:length(basis)) where T + return MomentMatrix(trimat(T, f, length(basis), σ), basis) end -MomentMatrix{T}(f::Function, x::AbstractVector, σ) where T = MomentMatrix{T}(f, monomial.(x), σ) -function MomentMatrix{T}(f::Function, x::AbstractVector) where T - σ, X = sortmonovec(x) - MomentMatrix{T}(f, X, σ) +function MomentMatrix{T}(f::Function, monos::AbstractVector) where T + σ, sorted_monos = sortmonovec(monos) + return MomentMatrix{T}(f, MB.MonomialBasis(sorted_monos), σ) end -MomentMatrix(f::Function, x) = MomentMatrix{Base.promote_op(f, Int, Int)}(f, x) -moment_matrix(f::Function, x) = MomentMatrix(f, x) """ moment_matrix(μ::Measure, x) @@ -51,27 +51,28 @@ function moment_matrix(μ::Measure{T}, X) where T end throw(ArgumentError("μ does not have the moment $(x)")) end - MomentMatrix{T}(getmom, X) + return MomentMatrix{T}(getmom, X) end -function MomentMatrix(Q::AbstractMatrix{T}, x, σ) where T - MomentMatrix{T}((i,j) -> Q[σ[i], σ[j]], x) +function MomentMatrix(Q::AbstractMatrix{T}, basis::MB.AbstractPolynomialBasis, σ) where T + return MomentMatrix{T}((i, j) -> Q[σ[i], σ[j]], basis) end -function MomentMatrix(Q::AbstractMatrix, x) - σ, X = sortmonovec(x) - MomentMatrix(Q, X, σ) +function MomentMatrix(Q::AbstractMatrix, monos::AbstractVector) + σ, sorted_monos = MP.sortmonovec(monos) + return MomentMatrix(Q, MB.MonomialBasis(sorted_monos), σ) end -moment_matrix(Q::AbstractMatrix, x) = MomentMatrix(Q, x) +moment_matrix(Q::AbstractMatrix, monos) = MomentMatrix(Q, monos) function getmat(μ::MomentMatrix) Matrix(μ.Q) end -function measure(ν::MomentMatrix) - n = length(ν.x) - measure(ν.Q.Q, [ν.x[i] * ν.x[j] for i in 1:n for j in 1:i]) +function measure(ν::MomentMatrix{T, <:MB.MonomialBasis}) where T + n = length(ν.basis) + monos = ν.basis.monomials + measure(ν.Q.Q, [monos[i] * monos[j] for i in 1:n for j in 1:i]) end -struct SparseMomentMatrix{T, MT <: MP.AbstractMonomial, MVT <: AbstractVector{MT}} <: AbstractMomentMatrix{T, MT, MVT} - sub_moment_matrices::Vector{MomentMatrix{T, MT, MVT}} +struct SparseMomentMatrix{T, B <: MB.AbstractPolynomialBasis} <: AbstractMomentMatrix{T, B} + sub_moment_matrices::Vector{MomentMatrix{T, B}} end diff --git a/test/moment_matrix.jl b/test/moment_matrix.jl index a8bfead..7e743f0 100644 --- a/test/moment_matrix.jl +++ b/test/moment_matrix.jl @@ -14,7 +14,7 @@ end μ = measure(1:3, [x^2, x*y, y^2]) X = [x, y] ν1 = moment_matrix(μ, X) - ν2 = moment_matrix((i, j) -> i + j - 1, X) + ν2 = MomentMatrix{Int}((i, j) -> i + j - 1, X) for ν in (ν1, ν2) @test ν.Q[1:4] == [1, 2, 2, 3] @test ν.Q[1, 1] == 1 @@ -27,7 +27,7 @@ end end @test_throws ArgumentError moment_matrix(measure([1], [x]), [y]) sparse_ν = SparseMomentMatrix([ν1, ν2]) - @test sparse_ν isa SparseMomentMatrix{Int,eltype(ν1.x),typeof(ν1.x)} + @test sparse_ν isa SparseMomentMatrix{Int,typeof(ν1.basis)} end # [HL05] Henrion, D. & Lasserre, J-B. @@ -65,7 +65,7 @@ end # β will be [z, y, x, y*z, x*z, x*y] x = monovec([z, y, x, z^2, y*z, y^2, x*z, x*y, x^2]) # x*β contains x*y*z, x^2*z, x^2*y which are not present so it show fail - V = MultivariateMoments.build_system(U', x, sqrt(eps(Float64))) + V = MultivariateMoments.build_system(U', MB.MonomialBasis(x), sqrt(eps(Float64))) @test iszerodimensional(V) testelements(V, [[2.0, 2.0, 2.0], [2.0, 2.0, 0.0], [2.0, 0.0, 2.0], [0.0, 2.0, 2.0], [2.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0], [0.0, 0.0, 0.0]], 1e-11) end diff --git a/test/runtests.jl b/test/runtests.jl index babc24b..4403f74 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,7 @@ -using MultivariateMoments using MultivariatePolynomials +import MultivariateBases +const MB = MultivariateBases +using MultivariateMoments using SemialgebraicSets using Test From 4e6195864fcdea507ec0b49facf19292210f2d7f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 27 Jan 2020 14:00:51 +0100 Subject: [PATCH 2/2] Add tests for map and copy for SymMatrix --- Project.toml | 1 + src/symmatrix.jl | 3 +++ test/moment_matrix.jl | 11 +++++++++-- 3 files changed, 13 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 07bfa7f..463bbd9 100644 --- a/Project.toml +++ b/Project.toml @@ -11,6 +11,7 @@ RowEchelon = "af85af4c-bcd5-5d23-b03a-a909639aa875" SemialgebraicSets = "8e049039-38e8-557d-ae3a-bc521ccf6204" [compat] +MultivariateBases = "0.1" MultivariatePolynomials = "0.3" RowEchelon = "0.1" SemialgebraicSets = "0.2" diff --git a/src/symmatrix.jl b/src/symmatrix.jl index cea13bf..0180a76 100644 --- a/src/symmatrix.jl +++ b/src/symmatrix.jl @@ -16,6 +16,9 @@ struct SymMatrix{T} <: AbstractMatrix{T} n::Int end +Base.copy(Q::SymMatrix) = SymMatrix(copy(Q.Q), Q.n) +Base.map(f::Function, Q::SymMatrix) = SymMatrix(map(f, Q.Q), Q.n) + # 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 diff --git a/test/moment_matrix.jl b/test/moment_matrix.jl index 7e743f0..157f81c 100644 --- a/test/moment_matrix.jl +++ b/test/moment_matrix.jl @@ -4,8 +4,15 @@ @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] + P = copy(Q) + @test P.n == 2 + symmetric_setindex!(P, 6, 2, 2) + @test Q.Q == [4, 5, 3] + @test P.n == 2 + @test P.Q == [4, 5, 6] + R = map(i -> i - 1, P) + @test R.n == 2 + @test R.Q == [3, 4, 5] end @testset "MomentMatrix" begin