Ortho and multi bases

src/MultivariateBases.jl

@@ -56,6 +56,8 @@ include("laguerre.jl")
 include("legendre.jl")
 include("chebyshev.jl")
 include("quotient.jl")
+include("orthonormal.jl")
+include("multi.jl")
 
 function algebra(
     basis::Union{QuotientBasis{Polynomial{B,M}},FullBasis{B,M},SubBasis{B,M}},

src/fixed.jl

@@ -1,66 +1,6 @@
-abstract type AbstractPolynomialVectorBasis{
-    PT<:MP.AbstractPolynomialLike,
-    PV<:AbstractVector{PT},
-} <: AbstractPolynomialBasis end
-
-function MP.monomial_type(
-    ::Type{<:AbstractPolynomialVectorBasis{PT}},
-) where {PT}
-    return MP.monomial_type(PT)
-end
-
-function empty_basis(
-    B::Type{<:AbstractPolynomialVectorBasis{PT,Vector{PT}}},
-) where {PT}
-    return B(PT[])
-end
-
-function MP.polynomial_type(
-    ::AbstractPolynomialVectorBasis{PT},
-    T::Type,
-) where {PT}
-    C = MP.coefficient_type(PT)
-    U = MA.promote_operation(*, C, T)
-    V = MA.promote_operation(+, U, U)
-    return MP.polynomial_type(PT, V)
-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{P},
-    ::Type{T},
-) where {P,T}
-    n = length(basis)
-    @assert size(Q) == (n, n)
-    PT = MP.polynomial_type(P, T)
-    return _convert(
-        PT,
-        mapreduce(
-            row -> begin
-                adjoint(basis.polynomials[row]) * mapreduce(
-                    col -> Q[row, col] * basis.polynomials[col],
-                    MA.add!!,
-                    1:n;
-                    init = MA.Zero(),
-                )
-            end,
-            MA.add!!,
-            1:n;
-            init = MA.Zero(),
-        ),
-    )::PT
-end
-
 """
-    struct FixedPolynomialBasis{PT<:MP.AbstractPolynomialLike, PV<:AbstractVector{PT}} <: AbstractPolynomialBasis
-        polynomials::PV
+    struct FixedBasis{T} <: SA.ExplicitBasis{T,I}
+        elements::Vector{T}
     end
 
 Polynomial basis with the polynomials of the vector `polynomials`.

src/multi.jl

@@ -0,0 +1,5 @@
+struct MultiBasis{T,I,B<:SA.ExplicitBasis{T,I}} <: SA.ExplicitBasis{T,I}
+    bases::Vector{B}
+end
+
+Base.length(basis::MultiBasis) = length(first(basis.bases))

src/orthonormal.jl

@@ -1,21 +1,34 @@
 """
-    struct OrthonormalCoefficientsBasis{PT<:MP.AbstractPolynomialLike, PV<:AbstractVector{PT}} <: AbstractPolynomialBasis
-        polynomials::PV
+    struct OrthonormalCoefficientsBasis{T,M<:AbstractMatrix{T},B} <: SA.ExplicitBasis
+        matrix::M
+        basis::B
     end
 
-Polynomial basis with the polynomials of the vector `polynomials` that are
-orthonormal with respect to the inner produce derived from the inner product
-of their coefficients.
-For instance, `FixedPolynomialBasis([1, x, 2x^2-1, 4x^3-3x])` is the Chebyshev
-polynomial basis for cubic polynomials in the variable `x`.
+Polynomial basis with the polynomials `algebra_element(matrix[i, :], basis)` for
+each row `i` of `matrix`. The basis is orthonormal with respect to the inner
+product derived from the inner product of their coefficients.
+For instance,
+```jldoctest
+julia> OrthonormalCoefficientsBasis(
+    [1 0 0 0
+     0 1 0 0
+    -1 0 2 0
+      -3   4],
+    SubBasis{Monomial}([1, x, x^2, x^3]),
+)
+```
+is the Chebyshev polynomial basis for cubic polynomials in the variable `x`.
 """
-struct OrthonormalCoefficientsBasis{
-    PT<:MP.AbstractPolynomialLike,
-    PV<:AbstractVector{PT},
-} <: AbstractBasis{PT,PV}
-    polynomials::PV
+struct OrthonormalCoefficientsBasis{T,B,M,V} <: SA.ExplicitBasis{
+    SA.AlgebraElement{Algebra{SubBasis{B,M,V},B,M},T,Vector{T}},
+    Int,
+}
+    matrix::Matrix{T}
+    basis::SubBasis{B,M,V}
 end
 
+Base.length(basis::OrthonormalCoefficientsBasis) = size(basis.matrix, 1)
+
 function LinearAlgebra.dot(
     p::MP.AbstractPolynomialLike{S},
     q::MP.AbstractPolynomialLike{T},

src/scaled.jl

@@ -51,7 +51,7 @@ end
 _float(::Type{T}) where {T<:Number} = float(T)
 # Could be for instance `MathOptInterface.ScalarAffineFunction{Float64}`
 # which does not implement `float`
-_float(::Type{F}) where {F} = F
+_float(::Type{F}) where {F} = MA.promote_operation(*, Float64, F)
 
 _promote_coef(::Type{T}, ::Type{ScaledMonomial}) where {T} = _float(T)
 

test/fixed.jl

@@ -1,9 +1,17 @@
 using Test
 using MultivariateBases
 using DynamicPolynomials
+import StarAlgebras as SA
 
 @polyvar x y
 
+struct ClassicalMul <: SA.MultiplicativeStructure end
+(::ClassicalMul)(a, b) = a * b
+
+gens = [algebra_element(1 - x^2), algebra_element(x^2 + 2)]
+basis = SA.FixedBasis(gens, ClassicalMul())
+a = SA.AlgebraElement([2, 3], basis)
+
 @testset "Polynomials" begin
     gens = [1 - x^2, x^2 + 2]
     basis = FixedPolynomialBasis(gens)