Lagrange basis (#47)

* Lagrange basis

* up

* up

* Fixes

* Fix format

* Add tests

* Fix format

* Add Random

* Fix doc

* Fixes

* Fix

Author: blegat

docs/src/index.md

@@ -11,7 +11,6 @@ end
 based on the [MultivariatePolynomials](https://github.com/JuliaAlgebra/MultivariatePolynomials.jl) API.
 
 ```@docs
-AbstractPolynomialBasis
 maxdegree_basis
 explicit_basis_covering
 ```

src/MultivariateBases.jl

@@ -4,7 +4,7 @@ import MutableArithmetics as MA
 import StarAlgebras as SA
 import MultivariatePolynomials as MP
 
-export AbstractPolynomialBasis, FullBasis, SubBasis
+export FullBasis, SubBasis
 export maxdegree_basis, explicit_basis_covering, empty_basis, monomial_index
 include("interface.jl")
 
@@ -55,15 +55,14 @@ export algebra_element,
     reccurence_second_coef,
     reccurence_third_coef,
     reccurence_deno_coef
-#include("fixed.jl")
 
 import LinearAlgebra
-#include("orthonormal.jl")
 include("orthogonal.jl")
 include("hermite.jl")
 include("laguerre.jl")
 include("legendre.jl")
 include("chebyshev.jl")
+include("lagrange.jl")
 include("quotient.jl")
 
 function algebra(

src/chebyshev.jl

@@ -89,6 +89,17 @@ function SA.coeffs!(
     return res
 end
 
+function transformation_to(
+    source::SubBasis{Chebyshev},
+    target::SubBasis{Monomial},
+)
+    A = zeros(Float64, length(target), length(source))
+    for (i, cheby) in enumerate(source)
+        A[:, i] = SA.coeffs(algebra_element(cheby), target)
+    end
+    return LinearAlgebra.UpperTriangular(A)
+end
+
 function SA.coeffs(
     cfs,
     source::SubBasis{Monomial},
@@ -97,17 +108,13 @@ function SA.coeffs(
     sub = explicit_basis_covering(target, source)
     # Need to make A square so that it's UpperTriangular
     extended = SubBasis{Monomial}(sub.monomials)
-    A = zeros(Float64, length(extended), length(sub))
-    for (i, cheby) in enumerate(sub)
-        A[:, i] = SA.coeffs(algebra_element(cheby), extended)
-    end
     ext = SA.coeffs(algebra_element(cfs, source), extended)
     return SA.SparseCoefficients(
         sub.monomials,
-        #LinearAlgebra.UpperTriangular(A) \ ext, # Julia v1.6 converts `A` to the eltype of the `result` which is bad for JuMP
+        #transformation_to(sub, extended) \ ext, # Julia v1.6 converts the matrix to the eltype of the `result` which is bad for JuMP
         LinearAlgebra.ldiv!(
-            zeros(_promote_coef(eltype(ext), Chebyshev), size(A, 2)),
-            LinearAlgebra.UpperTriangular(A),
+            zeros(_promote_coef(eltype(ext), Chebyshev), length(sub)),
+            transformation_to(sub, extended),
             ext,
         ),
     )
@@ -121,10 +128,7 @@ function SA.coeffs(cfs, ::FullBasis{Monomial}, target::FullBasis{Chebyshev})
     )
 end
 
-function degree_one_univariate_polynomial(
-    ::Type{Chebyshev},
-    variable::MP.AbstractVariable,
-)
+function degree_one_univariate_polynomial(::Type{Chebyshev}, variable)
     MA.@rewrite(variable + 0)
 end
 
@@ -146,10 +150,7 @@ Orthogonal polynomial with respect to the univariate weight function ``w(x) = \\
 """
 struct ChebyshevSecondKind <: AbstractChebyshev end
 
-function degree_one_univariate_polynomial(
-    ::Type{ChebyshevSecondKind},
-    variable::MP.AbstractVariable,
-)
+function degree_one_univariate_polynomial(::Type{ChebyshevSecondKind}, variable)
     MA.@rewrite(2variable + 0)
 end
 

src/fixed.jl

@@ -19,10 +19,7 @@ reccurence_first_coef(::Type{ProbabilistsHermite}, degree) = 1
 function reccurence_third_coef(::Type{ProbabilistsHermite}, degree)
     return -(degree - 1)
 end
-function degree_one_univariate_polynomial(
-    ::Type{ProbabilistsHermite},
-    variable::MP.AbstractVariable,
-)
+function degree_one_univariate_polynomial(::Type{ProbabilistsHermite}, variable)
     MA.@rewrite(1variable)
 end
 
@@ -47,10 +44,7 @@ Orthogonal polynomial with respect to the univariate weight function ``w(x) = \\
 struct PhysicistsHermite <: AbstractHermite end
 reccurence_first_coef(::Type{PhysicistsHermite}, degree) = 2
 reccurence_third_coef(::Type{PhysicistsHermite}, degree) = -2(degree - 1)
-function degree_one_univariate_polynomial(
-    ::Type{PhysicistsHermite},
-    variable::MP.AbstractVariable,
-)
+function degree_one_univariate_polynomial(::Type{PhysicistsHermite}, variable)
     MA.@rewrite(2variable)
 end
 

src/hermite.jl

@@ -1,24 +1,3 @@
-"""
-    abstract type AbstractPolynomialBasis end
-
-Polynomial basis of a subspace of the polynomials [Section~3.1.5, BPT12].
-
-[BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R.
-*Semidefinite Optimization and Convex Algebraic Geometry*.
-Society for Industrial and Applied Mathematics, **2012**.
-"""
-abstract type AbstractPolynomialBasis end
-
-# TODO breaking Should be underscore and only for internal use
-generators(basis::AbstractPolynomialBasis) = basis.polynomials
-
-function Base.getindex(
-    basis::AbstractPolynomialBasis,
-    I::AbstractVector{<:Integer},
-)
-    return typeof(basis)(generators(basis)[I])
-end
-
 """
     maxdegree_basis(basis::StarAlgebras.AbstractBasis, variables, maxdegree::Int)
 
@@ -27,7 +6,6 @@ Return the explicit version of `basis`generating all polynomials of degree up to
 """
 function maxdegree_basis end
 
-# TODO remove, not needed anymore
 """
     explicit_basis_covering(basis::StarAlgebras.AbstractBasis, target::StarAlgebras.ExplicitBasis)