From 0e595f45d78a348dec8d342be718a59231978bae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 20 Jun 2024 13:11:12 +0200 Subject: [PATCH] Ortho and multi bases --- src/MultivariateBases.jl | 2 ++ src/fixed.jl | 64 ++-------------------------------------- src/multi.jl | 5 ++++ src/orthonormal.jl | 37 +++++++++++++++-------- src/scaled.jl | 2 +- test/fixed.jl | 8 +++++ 6 files changed, 43 insertions(+), 75 deletions(-) create mode 100644 src/multi.jl diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl index 943268e..2d7334d 100644 --- a/src/MultivariateBases.jl +++ b/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}}, diff --git a/src/fixed.jl b/src/fixed.jl index 34ec967..070b81f 100644 --- a/src/fixed.jl +++ b/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`. diff --git a/src/multi.jl b/src/multi.jl new file mode 100644 index 0000000..ce1d1a4 --- /dev/null +++ b/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)) diff --git a/src/orthonormal.jl b/src/orthonormal.jl index 0827f76..fdcd441 100644 --- a/src/orthonormal.jl +++ b/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}, diff --git a/src/scaled.jl b/src/scaled.jl index 60403d9..4d8ec69 100644 --- a/src/scaled.jl +++ b/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) diff --git a/test/fixed.jl b/test/fixed.jl index d74cd3c..c090d80 100644 --- a/test/fixed.jl +++ b/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)