From 442cba8f705170b555088b7c12cb5a37ae1a08a0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 09:30:48 +0200 Subject: [PATCH 01/11] Lagrange basis --- src/lagrange.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 src/lagrange.jl diff --git a/src/lagrange.jl b/src/lagrange.jl new file mode 100644 index 0000000..744c61f --- /dev/null +++ b/src/lagrange.jl @@ -0,0 +1,10 @@ +struct ImplicitLagrangeBasis +end + +struct LagrangePolynomial{T,V<:AbstractVector{T}} + point::V +end + +struct LagrangeBasis{T,P<:AbstractVector{T},V<:AbstractVector{P}} + points::V +end From 2016b3f8900c42e816f6c56f19328d1e0ee07f60 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 11:12:04 +0200 Subject: [PATCH 02/11] up --- src/MultivariateBases.jl | 1 + src/chebyshev.jl | 21 +++++++++----- src/lagrange.jl | 61 ++++++++++++++++++++++++++++++++++++++-- src/orthogonal.jl | 15 ++++++++++ 4 files changed, 88 insertions(+), 10 deletions(-) diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl index c03ac5a..fcef0f7 100644 --- a/src/MultivariateBases.jl +++ b/src/MultivariateBases.jl @@ -64,6 +64,7 @@ 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..24aa106 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 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, ), ) diff --git a/src/lagrange.jl b/src/lagrange.jl index 744c61f..4788177 100644 --- a/src/lagrange.jl +++ b/src/lagrange.jl @@ -1,10 +1,65 @@ -struct ImplicitLagrangeBasis +# 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 + orthogonalize::Bool + function BoxSampling(lower, upper; sample_factor = 0, orthogonalize = false) + @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, orthogonalize) + end end -struct LagrangePolynomial{T,V<:AbstractVector{T}} +function sample(s::BoxSampling, n::Integer) where {T} + samples = rand(T, length(s.lower), n) .- inv(T(2)) + shift = (dom.u .+ dom.l) .* inv(T(2)) + for i in 1:n + for j in eachindex(s.lower) + samples[j, i] = samples[j, i] * (dom.u[j] - dom.l[j]) + shift[j] + end + end + return samples +end + +struct LagrangePolynomial{T,V} point::V end -struct LagrangeBasis{T,P<:AbstractVector{T},V<:AbstractVector{P}} +struct ImplicitLagrangeBasis{T,V,N<:AbstractNodes{T,V}} <: + SA.ImplicitBasis{LagrangePolynomial{T,V},V} + sampling::AbstractNodes{T,V} + function ImplicitLagrangeBasis(nodes::AbstractNodes{T,V}) where {T,V} + return new{T,V,typeof(nodes)}(nodes) + end +end +struct LagrangeBasis{T,P,V<:AbstractVector{P}} <: + SA.ExplicitBasis{LagrangePolynomial{T,V},Int} points::V end + +Base.length(basis::LagrangeBasis) = length(basis.points) + +function transformation_to(basis::SubBasis, lag::LagrangeBasis{T}) where {T} + V = Matrix{T}(undef, length(lag), length(basis)) + for +end + +function sample(s::AbstractNodes, basis::SubBasis) + samples = sample(s, length(basis) * s.sample_factor) + return eachcol(samples) +end + +function explicit_basis_covering( + implicit::ImplicitLagrangeBasis, + basis::SubBasis, +) + return LagrangeBasis(sample(implicit.sampling, length(basis))) +end diff --git a/src/orthogonal.jl b/src/orthogonal.jl index 9f5e0c6..093b416 100644 --- a/src/orthogonal.jl +++ b/src/orthogonal.jl @@ -62,6 +62,21 @@ 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 univariate_eval!( + values::Vector, + value, + ::Type{B}, +) where {B} + if 1 in eachindex(values) + values[1] = one(value) + end + if 2 in eachindex(values) + values[1] = one(value) + end + for d in 2:(length(values) - 1) + end +end + """ univariate_orthogonal_basis( B::Type{<:AbstractMultipleOrthogonal}, From f3197d68a83347678bcdc422662a8da673b1c86f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 12:10:53 +0200 Subject: [PATCH 03/11] up --- src/lagrange.jl | 73 ++++++++++++++++++++++++++++++++++++++++++----- src/monomial.jl | 6 ++++ src/orthogonal.jl | 54 ++++++++++++++++------------------- 3 files changed, 97 insertions(+), 36 deletions(-) diff --git a/src/lagrange.jl b/src/lagrange.jl index 4788177..22cc72d 100644 --- a/src/lagrange.jl +++ b/src/lagrange.jl @@ -18,12 +18,12 @@ mutable struct BoxSampling{T,V} <: AbstractNodes{T,V} end end -function sample(s::BoxSampling, n::Integer) where {T} +function sample(s::BoxSampling{T}, n::Integer) where {T} samples = rand(T, length(s.lower), n) .- inv(T(2)) - shift = (dom.u .+ dom.l) .* 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] * (dom.u[j] - dom.l[j]) + shift[j] + samples[j, i] = samples[j, i] * (s.upper[j] - s.lower[j]) + shift[j] end end return samples @@ -43,23 +43,82 @@ end struct LagrangeBasis{T,P,V<:AbstractVector{P}} <: SA.ExplicitBasis{LagrangePolynomial{T,V},Int} points::V + function LagrangeBasis(points::AbstractVector) + P = eltype(points) + return new{eltype(P), P, typeof(points)}(points) + end end Base.length(basis::LagrangeBasis) = length(basis.points) +function Base.getindex(basis::LagrangeBasis, I::AbstractVector{<:Integer}) + return LagrangeBasis(basis.points[I]) +end + +function eval_basis!(univariate_buffer, result, basis::SubBasis{B}, values) where {B} + for i in eachindex(values) + univariate_eval!(B, view(univariate_buffer, :, i), values[i]) + end + for i in eachindex(basis) + result[i] = one(eltype(result)) + exp = MP.exponents(basis.monomials[i]) + @assert length(exp) == length(values) + for j in eachindex(values) + result[i] = MA.operate!!(*, result[i], univariate_buffer[exp[j] + 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(basis)) V = Matrix{T}(undef, length(lag), length(basis)) - for + for i in eachindex(lag) + eval_basis!(univariate_buffer, view(V, i, :), basis, 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 function sample(s::AbstractNodes, basis::SubBasis) - samples = sample(s, length(basis) * s.sample_factor) - return eachcol(samples) + full = LagrangeBasis(eachcol(sample(s, num_samples(s.sample_factor, length(basis))))) + V = transformation_to(basis, full) + display(V) + F = LinearAlgebra.qr!(Matrix(V'), LinearAlgebra.ColumnNorm()) + display(F) + kept_indices = F.p[1:length(basis)] + return full[kept_indices] end function explicit_basis_covering( implicit::ImplicitLagrangeBasis, basis::SubBasis, ) - return LagrangeBasis(sample(implicit.sampling, length(basis))) + return sample(implicit.sampling, 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/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 093b416..cba17c4 100644 --- a/src/orthogonal.jl +++ b/src/orthogonal.jl @@ -62,19 +62,36 @@ 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 univariate_eval!( - values::Vector, +function recurrence_eval( + ::Type{B}, + previous::Vector, 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[1] = one(value) + 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 """ @@ -92,32 +109,11 @@ 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} From d5a88eef498f30c2624004777cb56a3cc19a45c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 14:22:25 +0200 Subject: [PATCH 04/11] Fixes --- src/lagrange.jl | 67 ++++++++++++++++++++++++++++++------------------- 1 file changed, 41 insertions(+), 26 deletions(-) diff --git a/src/lagrange.jl b/src/lagrange.jl index 22cc72d..d60b150 100644 --- a/src/lagrange.jl +++ b/src/lagrange.jl @@ -29,41 +29,57 @@ function sample(s::BoxSampling{T}, n::Integer) where {T} return samples end -struct LagrangePolynomial{T,V} - point::V +struct LagrangePolynomial{T,P,V} + variables::V + point::P end -struct ImplicitLagrangeBasis{T,V,N<:AbstractNodes{T,V}} <: - SA.ImplicitBasis{LagrangePolynomial{T,V},V} - sampling::AbstractNodes{T,V} - function ImplicitLagrangeBasis(nodes::AbstractNodes{T,V}) where {T,V} - return new{T,V,typeof(nodes)}(nodes) +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 -struct LagrangeBasis{T,P,V<:AbstractVector{P}} <: - SA.ExplicitBasis{LagrangePolynomial{T,V},Int} - points::V - function LagrangeBasis(points::AbstractVector) +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)}(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 Base.getindex(basis::LagrangeBasis, I::AbstractVector{<:Integer}) - return LagrangeBasis(basis.points[I]) + return LagrangeBasis(basis.variables, basis.points[I]) end -function eval_basis!(univariate_buffer, result, basis::SubBasis{B}, values) where {B} +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) - univariate_eval!(B, view(univariate_buffer, :, i), values[i]) + 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)) - exp = MP.exponents(basis.monomials[i]) - @assert length(exp) == length(values) for j in eachindex(values) - result[i] = MA.operate!!(*, result[i], univariate_buffer[exp[j] + 1, j]) + d = MP.degree(basis.monomials[i], variables[j]) + result[i] = MA.operate!!(*, result[i], univariate_buffer[d + 1, j]) end end return result @@ -72,10 +88,10 @@ 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(basis)) + 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, lag.points[i]) + eval_basis!(univariate_buffer, view(V, i, :), basis, MP.variables(lag), lag.points[i]) end return V end @@ -96,21 +112,20 @@ function num_samples(sample_factor, dim) return sample_factor * dim end -function sample(s::AbstractNodes, basis::SubBasis) - full = LagrangeBasis(eachcol(sample(s, num_samples(s.sample_factor, length(basis))))) +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) - display(V) F = LinearAlgebra.qr!(Matrix(V'), LinearAlgebra.ColumnNorm()) - display(F) kept_indices = F.p[1:length(basis)] - return full[kept_indices] + return LagrangeBasis(variables, eachcol(samples[:, kept_indices])) end function explicit_basis_covering( implicit::ImplicitLagrangeBasis, basis::SubBasis, ) - return sample(implicit.sampling, basis) + return sample(implicit.variables, implicit.nodes, basis) end function SA.coeffs(coeffs, source::SubBasis, target::LagrangeBasis) From e605773a4fa192b02fa36298025c63b4cbea6dc3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 14:23:49 +0200 Subject: [PATCH 05/11] Fix format --- src/lagrange.jl | 36 +++++++++++++++++++++++++++++------- src/orthogonal.jl | 17 +++++++++-------- 2 files changed, 38 insertions(+), 15 deletions(-) diff --git a/src/lagrange.jl b/src/lagrange.jl index d60b150..8afb2a5 100644 --- a/src/lagrange.jl +++ b/src/lagrange.jl @@ -38,7 +38,10 @@ 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} + function ImplicitLagrangeBasis( + variables, + nodes::AbstractNodes{T,P}, + ) where {T,P} return new{T,P,typeof(nodes),typeof(variables)}(variables, nodes) end end @@ -48,7 +51,10 @@ struct LagrangeBasis{T,P,U<:AbstractVector{P},V} <: points::U function LagrangeBasis(variables, points::AbstractVector) P = eltype(points) - return new{eltype(P),P,typeof(points),typeof(variables)}(variables, points) + return new{eltype(P),P,typeof(points),typeof(variables)}( + variables, + points, + ) end end @@ -59,16 +65,26 @@ function Base.getindex(basis::LagrangeBasis, I::AbstractVector{<:Integer}) return LagrangeBasis(basis.variables, basis.points[I]) end -function explicit_basis_type(::Type{<:ImplicitLagrangeBasis{T,_P,N,V}}) where {T,_P,N,V} +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} +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`") + 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) @@ -79,7 +95,7 @@ function eval_basis!(univariate_buffer, result, basis::SubBasis{B}, variables, v 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]) + result[i] = MA.operate!!(*, result[i], univariate_buffer[d+1, j]) end end return result @@ -91,7 +107,13 @@ function transformation_to(basis::SubBasis, lag::LagrangeBasis{T}) where {T} 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]) + eval_basis!( + univariate_buffer, + view(V, i, :), + basis, + MP.variables(lag), + lag.points[i], + ) end return V end diff --git a/src/orthogonal.jl b/src/orthogonal.jl index cba17c4..c9c88ef 100644 --- a/src/orthogonal.jl +++ b/src/orthogonal.jl @@ -77,19 +77,15 @@ function recurrence_eval( ) / d end -function univariate_eval!( - ::Type{B}, - values::AbstractVector, - value, -) where {B} +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) + for d in 2:(length(values)-1) + values[d+1] = recurrence_eval(B, view(values, 1:d), value, d) end return values end @@ -111,7 +107,12 @@ function univariate_orthogonal_basis( ) return univariate_eval!( B, - Vector{MP.polynomial_type(Polynomial{B,MP.monomial_type(variable)}, Int)}(undef, degree + 1), + Vector{ + MP.polynomial_type(Polynomial{B,MP.monomial_type(variable)}, Int), + }( + undef, + degree + 1, + ), variable, ) end From 44f9c61d604a5dd734c91824a1afe37724325dd4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 16:01:34 +0200 Subject: [PATCH 06/11] Add tests --- src/chebyshev.jl | 9 ++--- src/fixed.jl | 75 ----------------------------------------- src/hermite.jl | 4 +-- src/interface.jl | 22 ------------ src/lagrange.jl | 20 +++++++---- src/laguerre.jl | 2 +- src/legendre.jl | 2 +- src/orthogonal.jl | 2 +- src/orthonormal.jl | 48 -------------------------- test/fixed.jl | 80 ------------------------------------------- test/lagrange.jl | 61 +++++++++++++++++++++++++++++++++ test/orthonormal.jl | 82 --------------------------------------------- test/runtests.jl | 32 ++++-------------- 13 files changed, 89 insertions(+), 350 deletions(-) delete mode 100644 src/fixed.jl delete mode 100644 src/orthonormal.jl delete mode 100644 test/fixed.jl create mode 100644 test/lagrange.jl delete mode 100644 test/orthonormal.jl diff --git a/src/chebyshev.jl b/src/chebyshev.jl index 24aa106..21acd12 100644 --- a/src/chebyshev.jl +++ b/src/chebyshev.jl @@ -97,7 +97,7 @@ function transformation_to( for (i, cheby) in enumerate(source) A[:, i] = SA.coeffs(algebra_element(cheby), target) end - return A + return LinearAlgebra.UpperTriangular(A) end function SA.coeffs( @@ -128,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 @@ -155,7 +152,7 @@ struct ChebyshevSecondKind <: AbstractChebyshev end function degree_one_univariate_polynomial( ::Type{ChebyshevSecondKind}, - variable::MP.AbstractVariable, + 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..427df0d 100644 --- a/src/hermite.jl +++ b/src/hermite.jl @@ -21,7 +21,7 @@ function reccurence_third_coef(::Type{ProbabilistsHermite}, degree) end function degree_one_univariate_polynomial( ::Type{ProbabilistsHermite}, - variable::MP.AbstractVariable, + variable, ) MA.@rewrite(1variable) end @@ -49,7 +49,7 @@ 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, + 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 index 8afb2a5..f61e1d4 100644 --- a/src/lagrange.jl +++ b/src/lagrange.jl @@ -8,13 +8,12 @@ mutable struct BoxSampling{T,V} <: AbstractNodes{T,V} lower::V upper::V sample_factor::Int - orthogonalize::Bool - function BoxSampling(lower, upper; sample_factor = 0, orthogonalize = false) + 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, orthogonalize) + return new{eltype(V),V}(l, u, sample_factor) end end @@ -32,8 +31,12 @@ 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 @@ -45,6 +48,14 @@ struct ImplicitLagrangeBasis{T,P,N<:AbstractNodes{T,P},V} <: 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 @@ -61,9 +72,6 @@ end Base.length(basis::LagrangeBasis) = length(basis.points) MP.nvariables(basis::LagrangeBasis) = length(basis.variables) MP.variables(basis::LagrangeBasis) = basis.variables -function Base.getindex(basis::LagrangeBasis, I::AbstractVector{<:Integer}) - return LagrangeBasis(basis.variables, basis.points[I]) -end function explicit_basis_type( ::Type{<:ImplicitLagrangeBasis{T,_P,N,V}}, diff --git a/src/laguerre.jl b/src/laguerre.jl index c279409..d979f23 100644 --- a/src/laguerre.jl +++ b/src/laguerre.jl @@ -18,7 +18,7 @@ reccurence_deno_coef(::Type{Laguerre}, degree) = degree function degree_one_univariate_polynomial( ::Type{Laguerre}, - variable::MP.AbstractVariable, + variable, ) MA.@rewrite(1 - variable) end diff --git a/src/legendre.jl b/src/legendre.jl index 0595125..656a736 100644 --- a/src/legendre.jl +++ b/src/legendre.jl @@ -23,7 +23,7 @@ reccurence_deno_coef(::Type{Legendre}, degree) = degree function degree_one_univariate_polynomial( ::Type{Legendre}, - variable::MP.AbstractVariable, + variable, ) MA.@rewrite(variable + 0) end diff --git a/src/orthogonal.jl b/src/orthogonal.jl index c9c88ef..be12ca4 100644 --- a/src/orthogonal.jl +++ b/src/orthogonal.jl @@ -64,7 +64,7 @@ function reccurence_deno_coef end function recurrence_eval( ::Type{B}, - previous::Vector, + previous::AbstractVector, value, degree, ) where {B<:AbstractMultipleOrthogonal} 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/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..0911a5f --- /dev/null +++ b/test/lagrange.jl @@ -0,0 +1,61 @@ +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 From 23e4e989ade40af90046d7fc87a98a0ebc67b4ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 16:02:03 +0200 Subject: [PATCH 07/11] Fix format --- src/chebyshev.jl | 5 +---- src/hermite.jl | 10 ++-------- src/lagrange.jl | 15 +++++++++++---- src/laguerre.jl | 5 +---- src/legendre.jl | 5 +---- test/lagrange.jl | 19 +++++++++++-------- 6 files changed, 27 insertions(+), 32 deletions(-) diff --git a/src/chebyshev.jl b/src/chebyshev.jl index 21acd12..35619c1 100644 --- a/src/chebyshev.jl +++ b/src/chebyshev.jl @@ -150,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, -) +function degree_one_univariate_polynomial(::Type{ChebyshevSecondKind}, variable) MA.@rewrite(2variable + 0) end diff --git a/src/hermite.jl b/src/hermite.jl index 427df0d..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, -) +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, -) +function degree_one_univariate_polynomial(::Type{PhysicistsHermite}, variable) MA.@rewrite(2variable) end diff --git a/src/lagrange.jl b/src/lagrange.jl index f61e1d4..3547efd 100644 --- a/src/lagrange.jl +++ b/src/lagrange.jl @@ -32,11 +32,13 @@ struct LagrangePolynomial{T,P,V} variables::V point::P function LagrangePolynomial(variables, point) - return new{eltype(point),typeof(point),typeof(variables)}(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 @@ -49,9 +51,14 @@ struct ImplicitLagrangeBasis{T,P,N<:AbstractNodes{T,P},V} <: end end -function Base.getindex(implicit::ImplicitLagrangeBasis{T,P,N,V}, subs::Pair{V,P}) where {T,P,N,V} +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)`") + error( + "Variables `$(subs.first)` do not match Lagrange basis variables `$(implicit.variables)`", + ) end return LagrangePolynomial(implicit.variables, subs.second) end diff --git a/src/laguerre.jl b/src/laguerre.jl index d979f23..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, -) +function degree_one_univariate_polynomial(::Type{Laguerre}, variable) MA.@rewrite(1 - variable) end diff --git a/src/legendre.jl b/src/legendre.jl index 656a736..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, -) +function degree_one_univariate_polynomial(::Type{Legendre}, variable) MA.@rewrite(variable + 0) end diff --git a/test/lagrange.jl b/test/lagrange.jl index 0911a5f..fd17f08 100644 --- a/test/lagrange.jl +++ b/test/lagrange.jl @@ -8,23 +8,26 @@ 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]), - ) + implicit = + MB.ImplicitLagrangeBasis(x, MB.BoxSampling([-1, -1], UInt32[1, 1])) point = zeros(2) - poly = implicit[x => point] + 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]] + 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))}()) + 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]) From b338ab9f04941177a3163614a0a286af8a2c5ad4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 16:06:45 +0200 Subject: [PATCH 08/11] Add Random --- test/Project.toml | 1 + 1 file changed, 1 insertion(+) 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" From b69de57c91aef223a54413fea3dc022448a7b61f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 16:07:29 +0200 Subject: [PATCH 09/11] Fix doc --- docs/src/index.md | 1 - src/MultivariateBases.jl | 4 +--- 2 files changed, 1 insertion(+), 4 deletions(-) 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 fcef0f7..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,10 +55,8 @@ 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") From 6685c7b17b5125e7638650da714757519d1f3687 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 16:13:04 +0200 Subject: [PATCH 10/11] Fixes --- src/lagrange.jl | 13 ++++++++++--- test/lagrange.jl | 2 +- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/src/lagrange.jl b/src/lagrange.jl index 3547efd..51ba788 100644 --- a/src/lagrange.jl +++ b/src/lagrange.jl @@ -83,7 +83,7 @@ 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)) + points = _eachcol(ones(T, 1, 1)) P = eltype(points) return LagrangeBasis{eltype(P),P,typeof(points),V} end @@ -149,13 +149,20 @@ function num_samples(sample_factor, dim) return sample_factor * dim end +if VERSION >= v"1.10" + _eachcol(x) = eachcol(x) +else + # It is a `Base.Generator` so not an `AbstractVector` + _eachcol(x) = collect(eachcol(x)) +end + function sample(variables, s::AbstractNodes, basis::SubBasis) samples = sample(s, num_samples(s.sample_factor, length(basis))) - full = LagrangeBasis(variables, eachcol(samples)) + full = LagrangeBasis(variables, _eachcol(samples)) V = transformation_to(basis, full) F = LinearAlgebra.qr!(Matrix(V'), LinearAlgebra.ColumnNorm()) kept_indices = F.p[1:length(basis)] - return LagrangeBasis(variables, eachcol(samples[:, kept_indices])) + return LagrangeBasis(variables, _eachcol(samples[:, kept_indices])) end function explicit_basis_covering( diff --git a/test/lagrange.jl b/test/lagrange.jl index fd17f08..db988e7 100644 --- a/test/lagrange.jl +++ b/test/lagrange.jl @@ -28,7 +28,7 @@ function _test(B::Type) SA.SparseCoefficients(collect(monos), coeffs), MB.FullBasis{B,typeof(prod(x))}(), ) - @test SA.coeffs(a, lag) == SA.coeffs(coeffs, sub, lag) + @test SA.coeffs(a, lag) ≈ SA.coeffs(coeffs, sub, lag) @polyvar z bad = MB.SubBasis{B}([prod(x) * z]) err = ErrorException( From a77bc8076d1b77a7760f0ecb5b4e78da819962ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 Jul 2024 16:24:36 +0200 Subject: [PATCH 11/11] Fix --- src/lagrange.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/lagrange.jl b/src/lagrange.jl index 51ba788..7e4fcae 100644 --- a/src/lagrange.jl +++ b/src/lagrange.jl @@ -151,16 +151,18 @@ 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'), LinearAlgebra.ColumnNorm()) + F = LinearAlgebra.qr!(Matrix(V'), _column_norm()) kept_indices = F.p[1:length(basis)] return LagrangeBasis(variables, _eachcol(samples[:, kept_indices])) end