diff --git a/docs/src/index.md b/docs/src/index.md index e901ba5..9200511 100644 --- a/docs/src/index.md +++ b/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 ``` diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl index c03ac5a..8e35e58 100644 --- a/src/MultivariateBases.jl +++ b/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( diff --git a/src/chebyshev.jl b/src/chebyshev.jl index f94ab74..35619c1 100644 --- a/src/chebyshev.jl +++ b/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 diff --git a/src/fixed.jl b/src/fixed.jl deleted file mode 100644 index 34ec967..0000000 --- a/src/fixed.jl +++ /dev/null @@ -1,75 +0,0 @@ -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 - end - -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}, -} <: AbstractPolynomialVectorBasis{PT,PV} - polynomials::PV -end diff --git a/src/hermite.jl b/src/hermite.jl index a95353e..b98a041 100644 --- a/src/hermite.jl +++ b/src/hermite.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 diff --git a/src/interface.jl b/src/interface.jl index 2fe5c04..743a022 100644 --- a/src/interface.jl +++ b/src/interface.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) diff --git a/src/lagrange.jl b/src/lagrange.jl new file mode 100644 index 0000000..7e4fcae --- /dev/null +++ b/src/lagrange.jl @@ -0,0 +1,185 @@ +# Inspired by `Hypatia.jl/src/PolyUtils/realinterp.jl` + +export BoxSampling + +abstract type AbstractNodes{T,V} end + +mutable struct BoxSampling{T,V} <: AbstractNodes{T,V} + lower::V + upper::V + sample_factor::Int + function BoxSampling(lower, upper; sample_factor = 0) + @assert length(lower) == length(upper) + l = float.(lower) + u = float.(upper) + V = promote_type(typeof(l), typeof(u)) + return new{eltype(V),V}(l, u, sample_factor) + end +end + +function sample(s::BoxSampling{T}, n::Integer) where {T} + samples = rand(T, length(s.lower), n) .- inv(T(2)) + shift = (s.upper .+ s.lower) .* inv(T(2)) + for i in 1:n + for j in eachindex(s.lower) + samples[j, i] = samples[j, i] * (s.upper[j] - s.lower[j]) + shift[j] + end + end + return samples +end + +struct LagrangePolynomial{T,P,V} + variables::V + point::P + function LagrangePolynomial(variables, point) + return new{eltype(point),typeof(point),typeof(variables)}( + variables, + point, + ) + end +end + +struct ImplicitLagrangeBasis{T,P,N<:AbstractNodes{T,P},V} <: + SA.ImplicitBasis{LagrangePolynomial{T,P,V},Pair{V,P}} + variables::V + nodes::AbstractNodes{T,P} + function ImplicitLagrangeBasis( + variables, + nodes::AbstractNodes{T,P}, + ) where {T,P} + return new{T,P,typeof(nodes),typeof(variables)}(variables, nodes) + end +end + +function Base.getindex( + implicit::ImplicitLagrangeBasis{T,P,N,V}, + subs::Pair{V,P}, +) where {T,P,N,V} + if subs.first != implicit.variables + error( + "Variables `$(subs.first)` do not match Lagrange basis variables `$(implicit.variables)`", + ) + end + return LagrangePolynomial(implicit.variables, subs.second) +end + +struct LagrangeBasis{T,P,U<:AbstractVector{P},V} <: + SA.ExplicitBasis{LagrangePolynomial{T,P,V},Int} + variables::V + points::U + function LagrangeBasis(variables, points::AbstractVector) + P = eltype(points) + return new{eltype(P),P,typeof(points),typeof(variables)}( + variables, + points, + ) + end +end + +Base.length(basis::LagrangeBasis) = length(basis.points) +MP.nvariables(basis::LagrangeBasis) = length(basis.variables) +MP.variables(basis::LagrangeBasis) = basis.variables + +function explicit_basis_type( + ::Type{<:ImplicitLagrangeBasis{T,_P,N,V}}, +) where {T,_P,N,V} + points = _eachcol(ones(T, 1, 1)) + P = eltype(points) + return LagrangeBasis{eltype(P),P,typeof(points),V} +end + +function eval_basis!( + univariate_buffer, + result, + basis::SubBasis{B}, + variables, + values, +) where {B} + for v in MP.variables(basis) + if !(v in variables) + error( + "Cannot evaluate `$basis` as its variable `$v` is not part of the variables `$variables` of the `LagrangeBasis`", + ) + end + end + for i in eachindex(values) + l = MP.maxdegree(basis.monomials, variables[i]) + 1 + univariate_eval!(B, view(univariate_buffer, 1:l, i), values[i]) + end + for i in eachindex(basis) + result[i] = one(eltype(result)) + for j in eachindex(values) + d = MP.degree(basis.monomials[i], variables[j]) + result[i] = MA.operate!!(*, result[i], univariate_buffer[d+1, j]) + end + end + return result +end + +function transformation_to(basis::SubBasis, lag::LagrangeBasis{T}) where {T} + # To avoid allocating this too often, we allocate it once here + # and reuse it for each sample + univariate_buffer = Matrix{T}(undef, length(basis), MP.nvariables(lag)) + V = Matrix{T}(undef, length(lag), length(basis)) + for i in eachindex(lag) + eval_basis!( + univariate_buffer, + view(V, i, :), + basis, + MP.variables(lag), + lag.points[i], + ) + end + return V +end + +# Heuristic taken from Hypatia +function num_samples(sample_factor, dim) + if iszero(sample_factor) + if dim <= 12_000 + sample_factor = 10 + elseif dim <= 15_000 + sample_factor = 5 + elseif dim <= 22_000 + sample_factor = 2 + else + sample_factor = 1 + end + end + return sample_factor * dim +end + +if VERSION >= v"1.10" + _eachcol(x) = eachcol(x) + _column_norm() = LinearAlgebra.ColumnNorm() +else + # It is a `Base.Generator` so not an `AbstractVector` + _eachcol(x) = collect(eachcol(x)) + _column_norm() = Val(true) +end + +function sample(variables, s::AbstractNodes, basis::SubBasis) + samples = sample(s, num_samples(s.sample_factor, length(basis))) + full = LagrangeBasis(variables, _eachcol(samples)) + V = transformation_to(basis, full) + F = LinearAlgebra.qr!(Matrix(V'), _column_norm()) + kept_indices = F.p[1:length(basis)] + return LagrangeBasis(variables, _eachcol(samples[:, kept_indices])) +end + +function explicit_basis_covering( + implicit::ImplicitLagrangeBasis, + basis::SubBasis, +) + return sample(implicit.variables, implicit.nodes, basis) +end + +function SA.coeffs(coeffs, source::SubBasis, target::LagrangeBasis) + return transformation_to(source, target) * coeffs +end + +function SA.coeffs(coeffs, implicit::FullBasis, target::LagrangeBasis) + a = algebra_element(coeffs, implicit) + explicit = explicit_basis(a) + return SA.coeffs(SA.coeffs(a, explicit), explicit, target) +end diff --git a/src/laguerre.jl b/src/laguerre.jl index c279409..331cc20 100644 --- a/src/laguerre.jl +++ b/src/laguerre.jl @@ -16,10 +16,7 @@ reccurence_second_coef(::Type{Laguerre}, degree) = (2degree - 1) reccurence_third_coef(::Type{Laguerre}, degree) = -(degree - 1) reccurence_deno_coef(::Type{Laguerre}, degree) = degree -function degree_one_univariate_polynomial( - ::Type{Laguerre}, - variable::MP.AbstractVariable, -) +function degree_one_univariate_polynomial(::Type{Laguerre}, variable) MA.@rewrite(1 - variable) end diff --git a/src/legendre.jl b/src/legendre.jl index 0595125..3633de3 100644 --- a/src/legendre.jl +++ b/src/legendre.jl @@ -21,10 +21,7 @@ reccurence_first_coef(::Type{Legendre}, degree) = (2degree - 1) reccurence_third_coef(::Type{Legendre}, degree) = -(degree - 1) reccurence_deno_coef(::Type{Legendre}, degree) = degree -function degree_one_univariate_polynomial( - ::Type{Legendre}, - variable::MP.AbstractVariable, -) +function degree_one_univariate_polynomial(::Type{Legendre}, variable) MA.@rewrite(variable + 0) end diff --git a/src/monomial.jl b/src/monomial.jl index fd8d984..ffa1365 100644 --- a/src/monomial.jl +++ b/src/monomial.jl @@ -312,6 +312,12 @@ one get ths [`ScaledMonomial`](@ref). """ struct Monomial <: AbstractMonomial end +degree_one_univariate_polynomial(::Type{Monomial}, value) = value + +function recurrence_eval(::Type{Monomial}, previous, value, degree) + return previous[degree] * value +end + function (::Mul{Monomial})(a::MP.AbstractMonomial, b::MP.AbstractMonomial) return sparse_coefficients(a * b) end diff --git a/src/orthogonal.jl b/src/orthogonal.jl index 9f5e0c6..be12ca4 100644 --- a/src/orthogonal.jl +++ b/src/orthogonal.jl @@ -62,6 +62,34 @@ d_k p_k(x_i) = (a_k x_i + b_k) p_{k-1}(x_i) + c_k p_{k-2}(x_i) """ function reccurence_deno_coef end +function recurrence_eval( + ::Type{B}, + previous::AbstractVector, + value, + degree, +) where {B<:AbstractMultipleOrthogonal} + a = reccurence_first_coef(B, degree) + b = reccurence_second_coef(B, degree) + c = reccurence_third_coef(B, degree) + d = reccurence_deno_coef(B, degree) + return MA.@rewrite( + (a * value + b) * previous[degree] + c * previous[degree-1] + ) / d +end + +function univariate_eval!(::Type{B}, values::AbstractVector, value) where {B} + if 1 in eachindex(values) + values[1] = one(value) + end + if 2 in eachindex(values) + values[2] = degree_one_univariate_polynomial(B, value) + end + for d in 2:(length(values)-1) + values[d+1] = recurrence_eval(B, view(values, 1:d), value, d) + end + return values +end + """ univariate_orthogonal_basis( B::Type{<:AbstractMultipleOrthogonal}, @@ -77,32 +105,16 @@ function univariate_orthogonal_basis( variable::MP.AbstractVariable, degree::Integer, ) - @assert degree >= 0 - if degree == 0 - return MP.polynomial_type( - Polynomial{B,MP.monomial_type(variable)}, - Int, - )[one(variable)] - elseif degree == 1 - return push!( - univariate_orthogonal_basis(B, variable, 0), - degree_one_univariate_polynomial(B, variable), - ) - else - previous = univariate_orthogonal_basis(B, variable, degree - 1) - a = reccurence_first_coef(B, degree) - b = reccurence_second_coef(B, degree) - c = reccurence_third_coef(B, degree) - d = reccurence_deno_coef(B, degree) - next = MA.@rewrite( - (a * variable + b) * previous[degree] + c * previous[degree-1] - ) - if !isone(d) - next = next / d - end - push!(previous, next) - return previous - end + return univariate_eval!( + B, + Vector{ + MP.polynomial_type(Polynomial{B,MP.monomial_type(variable)}, Int), + }( + undef, + degree + 1, + ), + variable, + ) end function _covering(::FullBasis{B,M}, monos) where {B,M} diff --git a/src/orthonormal.jl b/src/orthonormal.jl deleted file mode 100644 index 0827f76..0000000 --- a/src/orthonormal.jl +++ /dev/null @@ -1,48 +0,0 @@ -""" - struct OrthonormalCoefficientsBasis{PT<:MP.AbstractPolynomialLike, PV<:AbstractVector{PT}} <: AbstractPolynomialBasis - polynomials::PV - 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`. -""" -struct OrthonormalCoefficientsBasis{ - PT<:MP.AbstractPolynomialLike, - PV<:AbstractVector{PT}, -} <: AbstractBasis{PT,PV} - polynomials::PV -end - -function LinearAlgebra.dot( - p::MP.AbstractPolynomialLike{S}, - q::MP.AbstractPolynomialLike{T}, - ::Type{<:OrthonormalCoefficientsBasis}, -) where {S,T} - s = zero(MA.promote_operation(*, S, T)) - terms_p = MP.terms(p) - terms_q = MP.terms(q) - tsp = iterate(terms_p) - tsq = iterate(terms_q) - while !isnothing(tsp) && !isnothing(tsq) - tp, sp = tsp - tq, sq = tsq - cmp = MP.compare(MP.monomial(tp), MP.monomial(tq)) - if iszero(cmp) - s += conj(MP.coefficient(tp)) * MP.coefficient(tq) - tsp = iterate(terms_p, sp) - tsq = iterate(terms_q, sq) - elseif cmp < 0 - tsp = iterate(terms_p, sp) - else - tsq = iterate(terms_q, sq) - end - end - return s -end - -function MP.coefficients(p, basis::OrthonormalCoefficientsBasis) - return [LinearAlgebra.dot(q, p, typeof(basis)) for q in generators(basis)] -end diff --git a/test/Project.toml b/test/Project.toml index 303677f..e253ba2 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -3,5 +3,6 @@ DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378" MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/fixed.jl b/test/fixed.jl deleted file mode 100644 index d74cd3c..0000000 --- a/test/fixed.jl +++ /dev/null @@ -1,80 +0,0 @@ -using Test -using MultivariateBases -using DynamicPolynomials - -@polyvar x y - -@testset "Polynomials" begin - 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 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]) - @test polynomial_type(basis, Int) == polynomial_type(x, Int) - @test polynomial(i -> i^2, basis) == 4x^2 + x -end -@testset "Terms" begin - basis = FixedPolynomialBasis([1, x, x^2]) - @test polynomial_type(basis, Int) == polynomial_type(x, Int) - @test polynomial(i -> (-1)^i, basis) == -x^2 + x - 1 -end -@testset "Linear" begin - basis = FixedPolynomialBasis([x, y]) - @test polynomial_type(basis, Int) == polynomial_type(x, Int) - @test polynomial(identity, basis) == x + 2y -end -@testset "One variable" begin - basis = FixedPolynomialBasis([x]) - @test polynomial_type(basis, Int) == polynomial_type(x, Int) - @test polynomial(i -> 5, basis) == 5x - @test typeof(polynomial(i -> 5, basis)) == polynomial_type(basis, Int) - @test typeof(polynomial(ones(Int, 1, 1), basis, Int)) <: - AbstractPolynomial{Int} - @test typeof(polynomial(ones(Int, 1, 1), basis, Float64)) <: - AbstractPolynomial{Float64} -end -@testset "Complex" begin - basis = FixedPolynomialBasis([(1 + 2im) * x]) - @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 - -@testset "Enumerate" begin - monos = [1, x, y^2] - basis = FixedPolynomialBasis(monos) - for (i, e) in enumerate(basis) - @test e == monos[i] - end -end diff --git a/test/lagrange.jl b/test/lagrange.jl new file mode 100644 index 0000000..db988e7 --- /dev/null +++ b/test/lagrange.jl @@ -0,0 +1,64 @@ +module TestLagrange + +using Random, Test +import StarAlgebras as SA +using DynamicPolynomials +import MultivariateBases as MB + +function _test(B::Type) + @polyvar x[1:2] + Random.seed!(0) + implicit = + MB.ImplicitLagrangeBasis(x, MB.BoxSampling([-1, -1], UInt32[1, 1])) + point = zeros(2) + poly = implicit[x=>point] + @test poly isa MB.LagrangePolynomial + @test poly.variables == x + @test poly.point === point + err = ErrorException( + "Variables `$([x[1]])` do not match Lagrange basis variables `$x`", + ) + @test_throws err implicit[[x[1]]=>[0.0]] + monos = monomials(x, 0:2) + coeffs = collect(eachindex(monos)) + sub = MB.SubBasis{B}(monos) + lag = MB.explicit_basis_covering(implicit, sub) + @test typeof(lag) == MB.explicit_basis_type(typeof(implicit)) + a = MB.algebra_element( + SA.SparseCoefficients(collect(monos), coeffs), + MB.FullBasis{B,typeof(prod(x))}(), + ) + @test SA.coeffs(a, lag) ≈ SA.coeffs(coeffs, sub, lag) + @polyvar z + bad = MB.SubBasis{B}([prod(x) * z]) + err = ErrorException( + "Cannot evaluate `$bad` as its variable `$z` is not part of the variables `$x` of the `LagrangeBasis`", + ) + @test_throws err SA.coeffs([1], bad, lag) +end + +test_mono() = _test(MB.Monomial) + +test_cheby() = _test(MB.Chebyshev) + +function test_num_samples() + @test MB.num_samples(1, 1) == 1 + @test MB.num_samples(0, 1) == 10 + @test MB.num_samples(0, 15_000) == 5 * 15_000 + @test MB.num_samples(0, 22_000) == 2 * 22_000 + @test MB.num_samples(0, 23_000) == 23_000 +end + +function runtests() + for name in names(@__MODULE__; all = true) + if startswith("$name", "test_") + @testset "$(name)" begin + getfield(@__MODULE__, name)() + end + end + end +end + +end + +TestLagrange.runtests() diff --git a/test/orthonormal.jl b/test/orthonormal.jl deleted file mode 100644 index 922ec0a..0000000 --- a/test/orthonormal.jl +++ /dev/null @@ -1,82 +0,0 @@ -using Test -using MultivariateBases -using DynamicPolynomials - -@polyvar x y - -@testset "Polynomials" begin - gens = [1 + x + y + x * y, 1 - x + y - x * y] / 2 - basis = OrthonormalCoefficientsBasis(gens) - @test iszero(dot(gens[1], gens[2], OrthonormalCoefficientsBasis)) - coefficient_test(basis, [2, -3]) - coefficient_test(basis, [-2im, 1 + 5im]) - coefficient_test(basis, [1im, 2im]) - @test polynomial_type(basis, Int) == polynomial_type(x, Float64) - @test polynomial(one, basis) == 1 + y - @test basis[1] == gens[1] - @test basis[2] == gens[2] - @test collect(basis) == gens - @test generators(basis) == gens - @test length(basis) == 2 - @test mindegree(basis) == 0 - @test mindegree(basis, x) == 0 - @test mindegree(basis, y) == 0 - @test maxdegree(basis) == 2 - @test maxdegree(basis, x) == 1 - @test maxdegree(basis, y) == 1 - @test extdegree(basis) == (0, 2) - @test extdegree(basis, x) == (0, 1) - @test extdegree(basis, y) == (0, 1) - @test variables(basis) == [x, y] - @test nvariables(basis) == 2 - @test sprint(show, basis) == - "OrthonormalCoefficientsBasis([0.5 + 0.5y + 0.5x + 0.5xy, 0.5 + 0.5y - 0.5x - 0.5xy])" - @test sprint(print, basis) == - "OrthonormalCoefficientsBasis([0.5 + 0.5*y + 0.5*x + 0.5*x*y, 0.5 + 0.5*y - 0.5*x - 0.5*x*y])" - b2 = basis[2:2] - @test length(b2) == 1 - @test b2[1] == gens[2] - b3 = basis[2:1] - @test isempty(b3) -end -@testset "Linear" begin - basis = OrthonormalCoefficientsBasis([x, y]) - @test polynomial_type(basis, Int) == polynomial_type(x, Int) - @test polynomial(identity, basis) == x + 2y -end -@testset "One variable" begin - basis = OrthonormalCoefficientsBasis([x]) - @test polynomial_type(basis, Int) == polynomial_type(x, Int) - @test polynomial(i -> 5, basis) == 5x - @test typeof(polynomial(i -> 5, basis)) == polynomial_type(basis, Int) - @test typeof(polynomial(ones(Int, 1, 1), basis, Int)) <: - AbstractPolynomial{Int} - @test typeof(polynomial(ones(Int, 1, 1), basis, Float64)) <: - AbstractPolynomial{Float64} -end -@testset "Complex" begin - for a in [1, -1, im, -im] - basis = OrthonormalCoefficientsBasis([a * x]) - @test 5x^2 == - @inferred polynomial(5ones(Int, 1, 1), basis, Complex{Int}) - @test 5x^2 == @inferred polynomial(5ones(Int, 1, 1), basis, Int) - coefficient_test(basis, [2]) - coefficient_test(basis, [-2im]) - coefficient_test(basis, [1 + 5im]) - end -end -@testset "Empty" begin - basis = OrthonormalCoefficientsBasis(typeof(x + 1)[]) - @test isempty(basis) - @test isempty(eachindex(basis)) - p = @inferred polynomial(zeros(Int, 0, 0), basis, Int) - @test iszero(p) -end - -@testset "Enumerate" begin - monos = [1, x, y^2] - basis = OrthonormalCoefficientsBasis(monos) - for (i, e) in enumerate(basis) - @test e == monos[i] - end -end diff --git a/test/runtests.jl b/test/runtests.jl index 38d0372..c276f31 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -240,30 +240,10 @@ function coefficient_test( return end -@testset "Monomial" begin - include("monomial.jl") -end -@testset "Scaled" begin - include("scaled.jl") -end -#@testset "Fixed" begin -# include("fixed.jl") -#end -#@testset "Orthonormal" begin -# include("orthonormal.jl") -#end -@testset "Hermite" begin - include("hermite.jl") -end -@testset "Laguerre" begin - include("laguerre.jl") -end -@testset "Legendre" begin - include("legendre.jl") -end -@testset "Chebyshev" begin - include("chebyshev.jl") -end -@testset "Quotient" begin - include("quotient.jl") +for file in readdir(@__DIR__) + if endswith(file, ".jl") && !in(file, ["runtests.jl"]) + @testset "$file" begin + include(joinpath(@__DIR__, file)) + end + end end