Add missing interface implementations for poly basis (#19)

* Add missing interface implementations for poly basis

* Fix format

* Bump MB to v0.5.3

* Fixes

* Add tests

Project.toml

@@ -8,13 +8,6 @@ MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 
 [compat]
-MultivariatePolynomials = "0.5"
+MultivariatePolynomials = "0.5.3"
 MutableArithmetics = "1"
 julia = "1.6"
-
-[extras]
-DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[targets]
-test = ["Test", "DynamicPolynomials"]

src/MultivariateBases.jl

@@ -22,7 +22,8 @@ export AbstractGegenbauerBasis,
     ChebyshevBasis,
     ChebyshevBasisFirstKind,
     ChebyshevBasisSecondKind
-export univariate_orthogonal_basis,
+export generators,
+    univariate_orthogonal_basis,
     reccurence_first_coef,
     reccurence_second_coef,
     reccurence_third_coef,

src/fixed.jl

@@ -3,23 +3,6 @@ abstract type AbstractPolynomialVectorBasis{
     PV<:AbstractVector{PT},
 } <: AbstractPolynomialBasis end
 
-Base.length(basis::AbstractPolynomialVectorBasis) = length(basis.polynomials)
-function Base.copy(basis::AbstractPolynomialVectorBasis)
-    return typeof(basis)(copy(basis.polynomials))
-end
-
-Base.firstindex(basis::AbstractPolynomialVectorBasis) = 1
-Base.lastindex(basis::AbstractPolynomialVectorBasis) = length(basis)
-function Base.getindex(basis::AbstractPolynomialVectorBasis, i::Int)
-    return basis.polynomials[i]
-end
-
-function MP.nvariables(basis::AbstractPolynomialVectorBasis)
-    return MP.nvariables(basis.polynomials)
-end
-function MP.variables(basis::AbstractPolynomialVectorBasis)
-    return MP.variables(basis.polynomials)
-end
 function MP.monomial_type(
     ::Type{<:AbstractPolynomialVectorBasis{PT}},
 ) where {PT}
@@ -40,36 +23,54 @@ function MP.polynomial_type(
     V = MA.promote_operation(+, U, U)
     return MP.polynomial_type(PT, V)
 end
-function MP.polynomial(f::Function, basis::AbstractPolynomialVectorBasis)
-    return MP.polynomial(
-        mapreduce(
-            ip -> f(ip[1]) * ip[2],
-            MA.add!!,
-            enumerate(basis.polynomials),
-        ),
-    )
+function MP.polynomial(
+    f::Function,
+    basis::AbstractPolynomialVectorBasis{P},
+) where {P}
+    if isempty(generators(basis))
+        return zero(P)
+    else
+        return MP.polynomial(
+            mapreduce(
+                ip -> f(ip[1]) * ip[2],
+                MA.add!!,
+                enumerate(basis.polynomials),
+            ),
+        )
+    end
 end
 
+function _poly(::MA.Zero, ::Type{P}, ::Type{T}) where {P,T}
+    return zero(MP.polynomial_type(P, T))
+end
+
+_convert(::Type{P}, p) where {P} = convert(P, p)
+_convert(::Type{P}, ::MA.Zero) where {P} = zero(P)
+
 function MP.polynomial(
     Q::AbstractMatrix,
-    basis::AbstractPolynomialVectorBasis,
-    T::Type,
-)
+    basis::AbstractPolynomialVectorBasis{P},
+    ::Type{T},
+) where {P,T}
     n = length(basis)
     @assert size(Q) == (n, n)
-    return MP.polynomial(
+    PT = MP.polynomial_type(P, T)
+    return _convert(
+        PT,
         mapreduce(
-            row ->
+            row -> begin
                 adjoint(basis.polynomials[row]) * mapreduce(
                     col -> Q[row, col] * basis.polynomials[col],
                     MA.add!!,
-                    1:n,
-                ),
+                    1:n;
+                    init = MA.Zero(),
+                )
+            end,
             MA.add!!,
-            1:n,
+            1:n;
+            init = MA.Zero(),
         ),
-        T,
-    )
+    )::PT
 end
 
 """

src/interface.jl

@@ -11,6 +11,42 @@ abstract type AbstractPolynomialBasis end
 
 generators(basis::AbstractPolynomialBasis) = basis.polynomials
 
+function Base.copy(basis::AbstractPolynomialBasis)
+    return typeof(basis)(copy(generators(basis)))
+end
+function Base.getindex(
+    basis::AbstractPolynomialBasis,
+    I::AbstractVector{<:Integer},
+)
+    return typeof(basis)(generators(basis)[I])
+end
+
+# Overload some of the `AbstractVector` interface for convenience
+Base.isempty(basis::AbstractPolynomialBasis) = isempty(generators(basis))
+Base.eachindex(basis::AbstractPolynomialBasis) = eachindex(generators(basis))
+Base.iterate(basis::AbstractPolynomialBasis) = iterate(generators(basis))
+Base.iterate(basis::AbstractPolynomialBasis, s) = iterate(generators(basis), s)
+Base.length(basis::AbstractPolynomialBasis) = length(generators(basis))
+Base.firstindex(basis::AbstractPolynomialBasis) = firstindex(generators(basis))
+Base.lastindex(basis::AbstractPolynomialBasis) = lastindex(generators(basis))
+Base.getindex(basis::AbstractPolynomialBasis, i::Int) = generators(basis)[i]
+
+# Overload some of the `MP` interface for convenience
+MP.mindegree(basis::AbstractPolynomialBasis) = MP.mindegree(generators(basis))
+MP.maxdegree(basis::AbstractPolynomialBasis) = MP.maxdegree(generators(basis))
+MP.extdegree(basis::AbstractPolynomialBasis) = MP.extdegree(generators(basis))
+function MP.mindegree(basis::AbstractPolynomialBasis, v)
+    return MP.mindegree(generators(basis), v)
+end
+function MP.maxdegree(basis::AbstractPolynomialBasis, v)
+    return MP.maxdegree(generators(basis), v)
+end
+function MP.extdegree(basis::AbstractPolynomialBasis, v)
+    return MP.extdegree(generators(basis), v)
+end
+MP.nvariables(basis::AbstractPolynomialBasis) = MP.nvariables(generators(basis))
+MP.variables(basis::AbstractPolynomialBasis) = MP.variables(generators(basis))
+
 function MP.polynomial(coefs::Vector, basis::AbstractPolynomialBasis)
     return MP.polynomial(i -> coefs[i], basis)
 end

src/monomial.jl

@@ -3,25 +3,22 @@ abstract type AbstractMonomialBasis{
     MV<:AbstractVector{MT},
 } <: AbstractPolynomialBasis end
 
-Base.length(basis::AbstractMonomialBasis) = length(basis.monomials)
-Base.copy(basis::AbstractMonomialBasis) = typeof(basis)(copy(basis.monomials))
-Base.firstindex(basis::AbstractMonomialBasis) = 1
-Base.lastindex(basis::AbstractMonomialBasis) = length(basis)
+generators(basis::AbstractMonomialBasis) = basis.monomials
 
-MP.nvariables(basis::AbstractMonomialBasis) = MP.nvariables(basis.monomials)
-MP.variables(basis::AbstractMonomialBasis) = MP.variables(basis.monomials)
 MP.monomial_type(::Type{<:AbstractMonomialBasis{MT}}) where {MT} = MT
 
 function empty_basis(MB::Type{<:AbstractMonomialBasis{MT}}) where {MT}
     return MB(MP.empty_monomial_vector(MT))
 end
+
 function maxdegree_basis(
     B::Type{<:AbstractMonomialBasis},
     variables,
     maxdegree::Int,
 )
     return B(MP.monomials(variables, 0:maxdegree))
 end
+
 function basis_covering_monomials(
     B::Type{<:AbstractMonomialBasis},
     monos::AbstractVector{<:MP.AbstractMonomial},
@@ -52,8 +49,6 @@ function merge_bases(basis1::MB, basis2::MB) where {MB<:AbstractMonomialBasis}
     return MB(monos), I1, I2
 end
 
-generators(basis::AbstractMonomialBasis) = basis.monomials
-
 """
     struct MonomialBasis{MT<:MP.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractPolynomialBasis
         monomials::MV
@@ -89,6 +84,7 @@ function MP.polynomial_type(
 ) where {MT}
     return MP.polynomial_type(MT, T)
 end
+
 MP.polynomial(f::Function, mb::MonomialBasis) = MP.polynomial(f, mb.monomials)
 
 function MP.polynomial(Q::AbstractMatrix, mb::MonomialBasis, T::Type)

test/Project.toml

@@ -0,0 +1,4 @@
+[deps]
+DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
+MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/fixed.jl

@@ -1,16 +1,36 @@
 using Test
 using MultivariateBases
 using DynamicPolynomials
+
 @polyvar x y
 
 @testset "Polynomials" begin
-    basis = FixedPolynomialBasis([1 - x^2, x^2 + 2])
+    gens = [1 - x^2, x^2 + 2]
+    basis = FixedPolynomialBasis(gens)
     @test polynomial_type(basis, Int) == polynomial_type(x, Int)
     @test polynomial(one, basis) == 3
     @test basis[1] == 1 - x^2
     @test basis[2] == x^2 + 2
+    @test collect(basis) == gens
+    @test generators(basis) == gens
+    @test length(basis) == 2
+    @test firstindex(basis) == 1
+    @test lastindex(basis) == 2
+    @test mindegree(basis) == 0
+    @test mindegree(basis, x) == 0
+    @test maxdegree(basis) == 2
+    @test maxdegree(basis, x) == 2
+    @test extdegree(basis) == (0, 2)
+    @test extdegree(basis, x) == (0, 2)
+    @test variables(basis) == [x]
+    @test nvariables(basis) == 1
     @test sprint(show, basis) == "FixedPolynomialBasis([1 - x², 2 + x²])"
     @test sprint(print, basis) == "FixedPolynomialBasis([1 - x^2, 2 + x^2])"
+    b2 = basis[2:2]
+    @test length(b2) == 1
+    @test b2[1] == x^2 + 2
+    b3 = basis[2:1]
+    @test isempty(b3)
 end
 @testset "Monomial" begin
     basis = FixedPolynomialBasis([x, x^2])
@@ -39,9 +59,16 @@ end
 end
 @testset "Complex" begin
     basis = FixedPolynomialBasis([(1 + 2im) * x])
-    @test 5x^2 == polynomial(ones(Int, 1, 1), basis, Complex{Int})
-    @test 5x^2 == polynomial(ones(Int, 1, 1), basis, Int)
+    @test 5x^2 == @inferred polynomial(ones(Int, 1, 1), basis, Complex{Int})
+    @test 5x^2 == @inferred polynomial(ones(Int, 1, 1), basis, Int)
     # TODO not inferred on Julia v1.0
     #@test 5x^2 == @inferred polynomial(ones(Int, 1, 1), basis, Complex{Int})
     #@test 5x^2 == @inferred polynomial(ones(Int, 1, 1), basis, Int)
 end
+@testset "Empty" begin
+    basis = FixedPolynomialBasis(typeof(x + 1)[])
+    @test isempty(basis)
+    @test isempty(eachindex(basis))
+    p = @inferred polynomial(zeros(Int, 0, 0), basis, Int)
+    @test iszero(p)
+end

test/monomial.jl

@@ -10,6 +10,13 @@ using DynamicPolynomials
     @test coefficients(x + 4y, MonomialBasis) == [4, 1]
     @test basis[1] == y
     @test basis[2] == x
+    @test generators(basis) == [y, x]
+    @test collect(basis) == [y, x]
+    @test length(basis) == 2
+    @test firstindex(basis) == 1
+    @test lastindex(basis) == 2
+    @test variables(basis) == [x, y]
+    @test nvariables(basis) == 2
     @test sprint(show, basis) == "MonomialBasis([y, x])"
     @test sprint(show, MIME"text/print"(), basis) == "MonomialBasis([y, x])"
     @test sprint(show, MIME"text/plain"(), basis) == "MonomialBasis([y, x])"
@@ -33,6 +40,15 @@ end
 @testset "merge_bases" begin
     basis1 = MonomialBasis([x^2, y^2])
     basis2 = MonomialBasis([x * y, y^2])
+    @test mindegree(basis2) == 2
+    @test mindegree(basis2, x) == 0
+    @test mindegree(basis2, y) == 1
+    @test maxdegree(basis2) == 2
+    @test maxdegree(basis2, x) == 1
+    @test maxdegree(basis2, y) == 2
+    @test extdegree(basis2) == (2, 2)
+    @test extdegree(basis2, x) == (0, 1)
+    @test extdegree(basis2, y) == (1, 2)
     basis, I1, I2 = MultivariateBases.merge_bases(basis1, basis2)
     @test basis.monomials == [y^2, x * y, x^2]
     @test I1 == [1, 0, 2]
@@ -42,3 +58,11 @@ end
 @testset "API degree = $degree" for degree in 0:3
     api_test(MonomialBasis, degree)
 end
+
+@testset "Empty" begin
+    basis = MonomialBasis(typeof(x^2)[])
+    @test isempty(basis)
+    @test isempty(eachindex(basis))
+    p = @inferred polynomial(zeros(Int, 0, 0), basis, Int)
+    @test iszero(p)
+end