Add Chebyshev basis

Project.toml

@@ -5,6 +5,12 @@ version = "0.1.0"
 
 [deps]
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
+MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
+
+[compat]
+MultivariatePolynomials = "0.3.5"
+MutableArithmetics = "0.2.1"
+julia = "1"
 
 [extras]
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"

src/MultivariateBases.jl

@@ -1,15 +1,21 @@
 module MultivariateBases
 
+import MutableArithmetics
+const MA = MutableArithmetics
+
 using MultivariatePolynomials
-const MP= MultivariatePolynomials
+const MP = MultivariatePolynomials
 
 export AbstractPolynomialBasis
-export FixedPolynomialBasis, ScaledMonomialBasis, MonomialBasis
-
+export maxdegree_basis, empty_basis
 include("interface.jl")
 
+export AbstractMonomialBasis, MonomialBasis, ScaledMonomialBasis
 include("monomial.jl")
-include("fixed.jl")
 include("scaled.jl")
 
+export FixedPolynomialBasis, ChebyshevBasis
+include("fixed.jl")
+include("chebyshev.jl")
+
 end # module

src/chebyshev.jl

@@ -0,0 +1,26 @@
+struct ChebyshevBasis{P} <: AbstractPolynomialVectorBasis{P, Vector{P}}
+    polynomials::Vector{P}
+end
+
+function chebyshev_polynomial_first_kind(variable::MP.AbstractVariable, degree::Integer)
+    @assert degree >= 0
+    if degree == 0
+        return [MA.@rewrite(0 * variable + 1)]
+    elseif degree == 1
+        return push!(chebyshev_polynomial_first_kind(variable, 0),
+                     MA.@rewrite(1 * variable + 0))
+    else
+        previous = chebyshev_polynomial_first_kind(variable, degree - 1)
+        next = MA.@rewrite(2variable * previous[degree] - previous[degree - 1])
+        push!(previous, next)
+        return previous
+    end
+end
+
+function maxdegree_basis(B::Type{ChebyshevBasis}, variables, maxdegree::Int)
+    univariate = [chebyshev_polynomial_first_kind(variable, maxdegree) for variable in variables]
+    return ChebyshevBasis([
+        prod(i -> univariate[i][degree(mono, variables[i]) + 1],
+             eachindex(variables))
+        for mono in MP.monomials(variables, 0:maxdegree)])
+end

src/fixed.jl

@@ -1,3 +1,33 @@
+abstract type AbstractPolynomialVectorBasis{PT<:MP.AbstractPolynomialLike, PV<:AbstractVector{PT}} <: AbstractPolynomialBasis end
+
+Base.length(basis::AbstractPolynomialVectorBasis) = length(basis.polynomials)
+Base.copy(basis::AbstractPolynomialVectorBasis) = typeof(basis)(copy(basis.polynomials))
+
+MP.nvariables(basis::AbstractPolynomialVectorBasis) = MP.nvariables(basis.polynomials)
+MP.variables(basis::AbstractPolynomialVectorBasis) = MP.variables(basis.polynomials)
+MP.monomialtype(::Type{<:AbstractPolynomialVectorBasis{PT}}) where PT = MP.monomialtype(PT)
+
+empty_basis(B::Type{<:AbstractPolynomialVectorBasis{PT, Vector{PT}}}) where PT = B(PT[])
+function MP.polynomialtype(basis::AbstractPolynomialVectorBasis{PT}, T::Type) where PT
+    C = MP.coefficienttype(PT)
+    U = MA.promote_operation(*, C, T)
+    V = MA.promote_operation(+, U, U)
+    return MP.polynomialtype(PT, V)
+end
+function MP.polynomial(f::Function, basis::AbstractPolynomialVectorBasis)
+    return mapreduce(ip -> f(ip[1]) * ip[2], MA.add!, enumerate(basis.polynomials))
+end
+
+function MP.polynomial(Q::AbstractMatrix, basis::AbstractPolynomialVectorBasis,
+                       T::Type)
+    n = length(basis)
+    @assert size(Q) == (n, n)
+    return mapreduce(row -> basis.polynomials[row] *
+        mapreduce(col -> Q[row, col] * basis.polynomials[col], MA.add!, 1:n),
+        MA.add!, 1:n)
+    return mapreduce()
+end
+
 """
     struct FixedPolynomialBasis{PT<:MP.AbstractPolynomialLike, PV<:AbstractVector{PT}} <: AbstractPolynomialBasis
         polynomials::PV
@@ -7,17 +37,6 @@ Polynomial basis with the polynomials of the vector `polynomials`.
 For instance, `FixedPolynomialBasis([1, x, 2x^2-1, 4x^3-3x])` is the Chebyshev
 polynomial basis for cubic polynomials in the variable `x`.
 """
-struct FixedPolynomialBasis{PT<:MP.AbstractPolynomialLike, PV<:AbstractVector{PT}} <: AbstractPolynomialBasis
+struct FixedPolynomialBasis{PT<:MP.AbstractPolynomialLike, PV<:AbstractVector{PT}} <: AbstractPolynomialVectorBasis{PT, PV}
     polynomials::PV
 end
-
-Base.length(basis::FixedPolynomialBasis) = length(basis.polynomials)
-empty_basis(::Type{FixedPolynomialBasis{PT, Vector{PT}}}) where PT = FixedPolynomialBasis(PT[])
-function MP.polynomialtype(mb::FixedPolynomialBasis{PT}, T::Type) where PT
-    C = MP.coefficienttype(PT)
-    U = typeof(zero(C) * zero(T) + zero(C) * zero(T))
-    MP.polynomialtype(PT, U)
-end
-function MP.polynomial(f::Function, fpb::FixedPolynomialBasis)
-    sum(ip -> f(ip[1]) * ip[2], enumerate(fpb.polynomials))
-end

src/monomial.jl

@@ -1,28 +1,15 @@
-"""
-    struct MonomialBasis{MT<:MP.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractPolynomialBasis
-        monomials::MV
-    end
+abstract type AbstractMonomialBasis{MT<:MP.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractPolynomialBasis end
 
-Monomial basis with the monomials of the vector `monomials`.
-For instance, `MonomialBasis([1, x, y, x^2, x*y, y^2])` is the monomial basis
-for the subspace of quadratic polynomials in the variables `x`, `y`.
-"""
-struct MonomialBasis{MT<:MP.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractPolynomialBasis
-    monomials::MV
-end
-MonomialBasis(monomials::AbstractVector) = MonomialBasis(MP.monovec(monomials))
-Base.copy(basis::MonomialBasis) = MonomialBasis(copy(basis.monomials))
+Base.length(basis::AbstractMonomialBasis) = length(basis.monomials)
+Base.copy(basis::AbstractMonomialBasis) = typeof(basis)(copy(basis.monomials))
 
-Base.length(basis::MonomialBasis) = length(basis.monomials)
-empty_basis(::Type{MonomialBasis{MT, MVT}}) where {MT, MVT} = MonomialBasis(MP.emptymonovec(MT))
-MP.nvariables(basis::MonomialBasis) = MP.nvariables(basis.monomials)
-MP.variables(basis::MonomialBasis) = MP.variables(basis.monomials)
-MP.monomialtype(::Type{<:MonomialBasis{MT}}) where MT = MT
-MP.polynomialtype(::Union{MonomialBasis{MT}, Type{<:MonomialBasis{MT}}}, T::Type) where MT = MP.polynomialtype(MT, T)
-MP.polynomial(f::Function, mb::MonomialBasis) = MP.polynomial(f, mb.monomials)
-MP.polynomial(Q::AbstractMatrix, mb::MonomialBasis, T::Type) = MP.polynomial(Q, mb.monomials, T)
-function MP.coefficients(p, ::Type{<:MonomialBasis})
-    return MP.coefficients(p)
+MP.nvariables(basis::AbstractMonomialBasis) = MP.nvariables(basis.monomials)
+MP.variables(basis::AbstractMonomialBasis) = MP.variables(basis.monomials)
+MP.monomialtype(::Type{<:AbstractMonomialBasis{MT}}) where MT = MT
+
+empty_basis(MB::Type{<:AbstractMonomialBasis{MT}}) where {MT} = MB(MP.emptymonovec(MT))
+function maxdegree_basis(B::Type{<:AbstractMonomialBasis}, variables, maxdegree::Int)
+    return B(MP.monomials(variables, 0:maxdegree))
 end
 
 # The `i`th index of output is the index of occurence of `x[i]` in `y`,
@@ -41,9 +28,32 @@ function multi_findsorted(x, y)
     return I
 end
 
-function merge_bases(basis1::MonomialBasis, basis2::MonomialBasis)
+function merge_bases(basis1::MB, basis2::MB) where MB<:AbstractMonomialBasis
     monos = MP.mergemonovec([basis1.monomials, basis2.monomials])
     I1 = multi_findsorted(monos, basis1.monomials)
     I2 = multi_findsorted(monos, basis2.monomials)
-    return monos, I1, I2
+    return MB(monos), I1, I2
+end
+
+"""
+    struct MonomialBasis{MT<:MP.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractPolynomialBasis
+        monomials::MV
+    end
+
+Monomial basis with the monomials of the vector `monomials`.
+For instance, `MonomialBasis([1, x, y, x^2, x*y, y^2])` is the monomial basis
+for the subspace of quadratic polynomials in the variables `x`, `y`.
+"""
+struct MonomialBasis{MT<:MP.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractMonomialBasis{MT, MV}
+    monomials::MV
+end
+MonomialBasis(monomials::AbstractVector) = MonomialBasis(MP.monovec(monomials))
+
+MP.polynomialtype(::Union{MonomialBasis{MT}, Type{<:MonomialBasis{MT}}}, T::Type) where MT = MP.polynomialtype(MT, T)
+MP.polynomial(f::Function, mb::MonomialBasis) = MP.polynomial(f, mb.monomials)
+
+MP.polynomial(Q::AbstractMatrix, mb::MonomialBasis, T::Type) = MP.polynomial(Q, mb.monomials, T)
+
+function MP.coefficients(p, ::Type{<:MonomialBasis})
+    return MP.coefficients(p)
 end

src/scaled.jl

@@ -20,16 +20,30 @@ Constraining the polynomial ``axy^2 + bxy`` to be zero with the scaled monomial
 *Semidefinite Optimization and Convex Algebraic Geometry*.
 Society for Industrial and Applied Mathematics, **2012**.
 """
-struct ScaledMonomialBasis{MT<:MP.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractPolynomialBasis
+struct ScaledMonomialBasis{MT<:MP.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractMonomialBasis{MT, MV}
     monomials::MV
 end
 ScaledMonomialBasis(monomials) = ScaledMonomialBasis(monovec(monomials))
 
-Base.length(basis::ScaledMonomialBasis) = length(basis.monomials)
-empty_basis(::Type{ScaledMonomialBasis{MT, MVT}}) where {MT, MVT} = ScaledMonomialBasis(MP.emptymonovec(MT))
-MP.polynomialtype(mb::ScaledMonomialBasis{MT}, T::Type) where MT = MP.polynomialtype(MT, promote_type(T, Float64))
+MP.polynomialtype(::ScaledMonomialBasis{MT}, T::Type) where MT = MP.polynomialtype(MT, promote_type(T, Float64))
+MP.polynomial(f::Function, basis::ScaledMonomialBasis) = MP.polynomial(i -> scaling(basis.monomials[i]) * f(i), basis.monomials)
+
+function Base.promote_rule(::Type{ScaledMonomialBasis{MT, MV}}, ::Type{MonomialBasis{MT, MV}}) where {MT, MV}
+    return MonomialBasis{MT, MV}
+end
+
+function change_basis(Q::AbstractMatrix, basis::ScaledMonomialBasis{MT, MV}, B::Type{MonomialBasis{MT, MV}}) where {MT, MV}
+    n = length(basis)
+    scalings = map(scaling, basis.monomials)
+    scaled_Q = [Q[i, j] * scalings[i] * scalings[j] for i in 1:n, j in 1:n]
+    return scaled_Q, MonomialBasis(basis.monomials)
+end
+
+function MP.polynomial(Q::AbstractMatrix, basis::ScaledMonomialBasis{MT, MV}, T::Type) where {MT, MV}
+    return MP.polynomial(change_basis(Q, basis, MonomialBasis{MT, MV})..., T)
+end
+
 scaling(m::MP.AbstractMonomial) = √(factorial(MP.degree(m)) / prod(factorial, MP.exponents(m)))
-MP.polynomial(f::Function, mb::ScaledMonomialBasis) = MP.polynomial(i -> scaling(mb.monomials[i]) * f(i), mb.monomials)
 unscale_coef(t::MP.AbstractTerm) = coefficient(t) / scaling(monomial(t))
 function MP.coefficients(p, ::Type{<:ScaledMonomialBasis})
     return unscale_coef.(MP.terms(p))

test/chebyshev.jl

@@ -0,0 +1,16 @@
+using Test
+using MultivariateBases
+using DynamicPolynomials
+@polyvar x y
+
+@testset "Univariate" begin
+    basis = maxdegree_basis(ChebyshevBasis, (x,), 3)
+    @test basis.polynomials[4] == 1
+    @test basis.polynomials[3] == x
+    @test basis.polynomials[2] == 2x^2 - 1
+    @test basis.polynomials[1] == 4x^3 - 3x
+end
+
+@testset "API degree = $degree" for degree in 0:3
+    api_test(ChebyshevBasis, degree)
+end

test/monomial.jl

@@ -22,3 +22,6 @@ end
     @test polynomial(i -> i^2, basis) == x^2 + 4x*y + 9y^2
     @test coefficients(x^2 + 4x*y + 9y^2, MonomialBasis) == [1, 4, 9]
 end
+@testset "API degree = $degree" for degree in 0:3
+    api_test(MonomialBasis, degree)
+end

test/runtests.jl

@@ -2,12 +2,34 @@ using Test
 
 using MultivariateBases
 
+using DynamicPolynomials
+
+function api_test(B::Type{<:AbstractPolynomialBasis}, degree)
+    @polyvar x[1:2]
+    basis = maxdegree_basis(B, x, degree)
+    n = binomial(2 + degree, 2)
+    @test length(basis) == n
+    @test typeof(copy(basis)) == typeof(basis)
+    @test nvariables(basis) == 2
+    @test variables(basis) == x
+    @test monomialtype(typeof(basis)) == monomialtype(x[1])
+    @test typeof(empty_basis(typeof(basis))) == typeof(basis)
+    @test length(empty_basis(typeof(basis))) == 0
+    @test polynomialtype(basis, Float64) == polynomialtype(x[1], Float64)
+    @test polynomial(i -> 0.0, basis) isa polynomialtype(basis, Float64)
+    @test polynomial(zeros(n, n), basis, Float64) isa polynomialtype(basis, Float64)
+    @test polynomial(ones(n, n), basis, Float64) isa polynomialtype(basis, Float64)
+end
+
 @testset "Monomial" begin
     include("monomial.jl")
 end
+@testset "Scaled" begin
+    include("scaled.jl")
+end
 @testset "Fixed" begin
     include("fixed.jl")
 end
-@testset "Scaled" begin
-    include("scaled.jl")
+@testset "Chebyshev" begin
+    include("chebyshev.jl")
 end

test/scaled.jl

@@ -22,3 +22,6 @@ end
     @test polynomial(i -> i^2, basis) == x^2 + 4*√2*x*y + 9y^2
     @test coefficients(x^2 + 4x*y + 9y^2, ScaledMonomialBasis) == [1, 4 / √2, 9]
 end
+@testset "API degree = $degree" for degree in 0:3
+    api_test(ScaledMonomialBasis, degree)
+end