Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

working with tw/MB #16

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/atomic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -45,9 +45,9 @@ end

measure(η::AtomicMeasure, x) = Measure(η, x)
function Measure(η::AtomicMeasure{T}, x::AbstractVector{TT}) where {T, TT}
Measure{T, monomialtype(TT), monovectype(x)}(η, x)
Measure{T, MB.MonomialBasis}(η, x)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The full type is MB.MonomialBasis{...}

end
function Measure{T, MT, MVT}(η::AtomicMeasure{T}, x) where {T, MT, MVT}
function Measure{T, BT}(η::AtomicMeasure{T}, x) where {T, BT}
X = monovec(x)
sum(atom.weight * dirac(X, η.variables => atom.center) for atom in η.atoms)
end
Expand Down
10 changes: 5 additions & 5 deletions src/expectation.jl
Original file line number Diff line number Diff line change
@@ -1,16 +1,16 @@
const APL = AbstractPolynomialLike

function _expectation(μ::Measure, p::APL, f)
function _expectation(μ::Measure{T, BT}, p::APL, f) where {T, BT}
i = 1
s = 0
for t in terms(p)
while i <= length(μ.x) && monomial(t) != μ.x[i]
for (c, m) in zip(coefficients(p, BT), MB.basis_covering_monomials(BT, monomials(p)))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

coefficients(p, BT) computes maxdegree_basis(Basis, variables(p), maxdegree(p)) which may not match with the basis basis = MB.basis_covering_monomials(BT, monomials(p)), use coefficients(p, basis) instead maybe

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See my comment in MultivariateBases.

while i <= length(μ.x) && m != μ.x[i]
i += 1
end
if i > length(μ.x)
error("The polynomial $p has a nonzero term $t with monomial $(t.x) for which the expectation is not known in $μ")
error("The polynomial $p has a nonzero term $(c*m) with basis function $(m) for which the expectation is not known in $μ")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would not call it term $(c * m) but rather coordinate $c or something like that

end
s += f(μ.a[i], coefficient(t))
s += f(μ.a[i], c)
i += 1
end
s
Expand Down
4 changes: 2 additions & 2 deletions src/extract.jl
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ function build_system(U::AbstractMatrix, basis::MB.MonomialBasis, ztol)
r = length(pivots)
U = U[keep, :]
end
monos = basis.monomials
monos = basis.elements
β = 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
Expand Down Expand Up @@ -181,7 +181,7 @@ function extractatoms(ν::MomentMatrix{T}, ranktol, args...) where T
vars = variables(μ)
A = Matrix{T}(undef, length(μ.x), r)
for i in 1:r
A[:, i] = dirac(μ.x, vars => centers[i]).a
A[:, i] = dirac(μ.x.elements, vars => centers[i]).a
end
weights = A \ μ.a
isf = isfinite.(weights)
Expand Down
27 changes: 20 additions & 7 deletions src/measure.jl
Original file line number Diff line number Diff line change
@@ -1,17 +1,22 @@
export measure, dirac
export variables, monomials, moments
export variables, base_functions, moments

# If a monomial is not in x, it does not mean that the moment is zero, it means that it is unknown/undefined
struct Measure{T, MT <: AbstractMonomial, MVT <: AbstractVector{MT}} <: AbstractMeasure{T}
struct Measure{T, B<:MB.AbstractPolynomialBasis} <: AbstractMeasure{T}
a::Vector{T}
x::MVT

function Measure{T, MT, MVT}(a::Vector{T}, x::MVT) where {T, MT, MVT}
x::B
blegat marked this conversation as resolved.
Show resolved Hide resolved
function Measure(a::Vector{T}, x::B) where {T, B<:MB.AbstractPolynomialBasis}
@assert length(a) == length(x)
new(a, x)
new{T, B}(a, x)
end
end
Measure(a::AbstractVector{T}, x::AbstractVector{TT}) where {T, TT <: AbstractTermLike} = Measure{T, monomialtype(TT), monovectype(x)}(monovec(a, x)...)

basis_functions_type(::Measure{T, BT}) where {T,BT} = BT

function Measure(a::AbstractVector{T}, x::AbstractVector{TT}) where {T, TT <: AbstractTermLike}
b, y = monovec(a, x)
return Measure(b, MB.MonomialBasis(y))
end

"""
measure(a, X::AbstractVector{<:AbstractMonomial})
Expand All @@ -34,6 +39,14 @@ Returns an iterator over the monomials of `μ` sorted in the decreasing order.
"""
MP.monomials(μ::Measure) = μ.x

"""
base_functions(μ::AbstractMeasureLike)

Returns an iterator over the base_functions of `μ` sorted in the decreasing order.
"""
base_functions(μ::Measure) = μ.x


"""
moments(μ::AbstractMeasureLike)

Expand Down
29 changes: 21 additions & 8 deletions src/moment.jl
Original file line number Diff line number Diff line change
@@ -1,16 +1,16 @@
export moment, moment_value, monomial
export moment, moment_value, monomial, base_function

struct Moment{T, MT <: AbstractMonomial} <: AbstractMoment{T}
struct Moment{T, PT <: AbstractPolynomialLike} <: AbstractMoment{T}
α::T
x::MT
x::PT
end

"""
moment(α, m::AbstractMonomial)
moment(α, m::AbstractPolynomialLike)

Creates the moment of the monomial `m` of value `α`.
Creates the moment of the base function `m` of value `α`.
"""
moment(α, m::AbstractMonomial) = Moment(α, m)
moment(α, m::AbstractPolynomialLike) = Moment(α, m)

"""
moment_value(m::AbstractMomentLike)
Expand All @@ -34,8 +34,21 @@ Calling `monomial(moment(3.1, x*y^2))` should return `x*y^2`.
"""
MP.monomial(m::Moment) = m.x

for f in [:variables, :nvariables, :exponents, :degree, :powers]
"""
base_function(m::AbstractMomentLike)

Returns the base_function of the moment `m`.

## Examples

Calling `base_function(moment(3.1, x*y^2))` should return `x*y^2`.
"""
base_function(m::Moment) = m.x

for f in [:variables, :nvariables, :exponents, :powers]
@eval begin
MP.$f(m::AbstractMomentLike) = $f(monomial(m))
MP.$f(m::AbstractMomentLike) = $f(base_function(m))
end
end

MP.degree(m::AbstractMomentLike) = MP.maxdegree(base_function(m))
2 changes: 1 addition & 1 deletion src/moment_matrix.jl
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ end

function measure(ν::MomentMatrix{T, <:MB.MonomialBasis}) where T
n = length(ν.basis)
monos = ν.basis.monomials
monos = ν.basis.elements
measure(ν.Q.Q, [monos[i] * monos[j] for i in 1:n for j in 1:i])
end

Expand Down
11 changes: 11 additions & 0 deletions test/expectation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -7,4 +7,15 @@
@test_throws ErrorException dot(x[1] * x[2] * x[3], m)
@test dot(0.5 * x[1] * x[2]^2, m) == 2.0
@test dot(m, x[1] * x[3]) == 3

@testset "Expectation - MB" begin
Mod.@polyvar x[1:2]
p = x[2] - 2x[1]*x[2]^2 + 3x[2]*x[1] - 5x[1]^3
m = measure([1, 0, 2, 3, 0, 4, 0, 2, 5, 0], maxdegree_basis(ChebyshevBasis, x, 3))
@test MultivariateMoments.expectation(m, p) == MultivariateMoments.expectation(p, m) == 4.25
@test_throws ErrorException dot((x[1] * x[2])^2, m)
@test dot(0.5 * x[1] * x[2]^2, m) == 1.0
@test dot(m, x[1] * x[2]) == 4
end

end
21 changes: 19 additions & 2 deletions test/measure.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,8 @@
@test_throws ArgumentError measure([1, 2], [x, x*y, y])
@test_throws ArgumentError measure([1, 2, 3, 4], [x, x*y, y])
μ = measure([1, 0, 2, 3], [x^2*y^2, y*x^2, x*y*x^2, x*y^2])
@test monomials(μ) == [x^3*y, x^2*y^2, x^2*y, x*y^2]
@test monomial.(moments(μ)) == [x^3*y, x^2*y^2, x^2*y, x*y^2]
@test base_functions(μ) == MonomialBasis(monovec([x^3*y, x^2*y^2, x^2*y, x*y^2]))
@test base_function.(moments(μ)) == [x^3*y, x^2*y^2, x^2*y, x*y^2]
@test MultivariateMoments.moment_value.(moments(μ)) == [2, 1, 0, 3]
@test all(nvariables.(moments(μ)) .== 2)
@test degree.(moments(μ)) == [4, 4, 3, 3]
Expand All @@ -13,4 +13,21 @@
@test (2 * μ).a == [4, 2, 0, 6]
@test (μ * 3).a == [6, 3, 0, 9]
#@test_throws ArgumentError measure([1], [x]) + measure([1], [y])

@testset "Measure MB" begin
Mod.@polyvar x y
@test_throws AssertionError measure([1, 0, 2, 3, 0], maxdegree_basis(ChebyshevBasis, [x, y], 2))
@test_throws AssertionError measure([1, 0, 2, 3, 0, 4], maxdegree_basis(ChebyshevBasis, [x, y], 3))
μ = measure([1, 0, 2, 3, 0, 4], maxdegree_basis(ChebyshevBasis, [x, y], 2))
@test base_functions(μ) == maxdegree_basis(ChebyshevBasis, [x, y], 2)
@test base_function.(moments(μ)) == [2.0x^2 - 1.0, x*y, 2.0y^2 - 1.0, x, y, 1.0]
@test MultivariateMoments.moment_value.(moments(μ)) == [1, 0, 2, 3, 0, 4]
@test all(nvariables.(moments(μ)) .== 2)
@test degree.(moments(μ)) == [2, 2, 2, 1, 1, 0]
@test μ.a == [1, 0, 2, 3, 0, 4]
@test (-μ).a == -[1, 0, 2, 3, 0, 4]
@test (2 * μ).a == [1, 0, 2, 3, 0, 4].*2
@test (μ * 3).a == [1, 0, 2, 3, 0, 4].*3

end
end
2 changes: 1 addition & 1 deletion test/runtests.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
using MultivariatePolynomials
import MultivariateBases
using MultivariateBases
const MB = MultivariateBases
using MultivariateMoments
using SemialgebraicSets
Expand Down