From fee4ab8c212b4818ecad9f1a0e5dc33a61834566 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 14 Jun 2024 17:23:51 +0200 Subject: [PATCH] Fix format --- src/MultivariateBases.jl | 15 +++++----- src/chebyshev.jl | 38 +++++++++++++++++++------ src/monomial.jl | 60 ++++++++++++++++++++++++++++------------ src/orthogonal.jl | 20 ++++++++++---- src/polynomial.jl | 32 ++++++++++++++------- src/quotient.jl | 37 ++++++++++++++++++++----- src/scaled.jl | 45 ++++++++++++++++++++++++------ test/chebyshev.jl | 3 +- 8 files changed, 186 insertions(+), 64 deletions(-) diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl index b468245..bea4086 100644 --- a/src/MultivariateBases.jl +++ b/src/MultivariateBases.jl @@ -14,7 +14,8 @@ MP.monomial_type(::Type{<:SA.AlgebraElement{A}}) where {A} = MP.monomial_type(A) function MP.polynomial_type(::Type{<:SA.AlgebraElement{A,T}}) where {A,T} return MP.polynomial_type(A, T) end -struct Algebra{BT,B,M} <: SA.AbstractStarAlgebra{Polynomial{B,M},Polynomial{B,M}} +struct Algebra{BT,B,M} <: + SA.AbstractStarAlgebra{Polynomial{B,M},Polynomial{B,M}} basis::BT end MP.monomial_type(::Type{<:Algebra{B}}) where {B} = MP.monomial_type(B) @@ -56,17 +57,17 @@ include("legendre.jl") include("chebyshev.jl") include("quotient.jl") -function algebra(basis::Union{QuotientBasis{Polynomial{B,M}},FullBasis{B,M},SubBasis{B,M}}) where {B,M} +function algebra( + basis::Union{QuotientBasis{Polynomial{B,M}},FullBasis{B,M},SubBasis{B,M}}, +) where {B,M} return Algebra{typeof(basis),B,M}(basis) end function MA.promote_operation( ::typeof(algebra), - BT::Type{<:Union{ - QuotientBasis{Polynomial{B,M}}, - FullBasis{B,M}, - SubBasis{B,M}, - }}, + BT::Type{ + <:Union{QuotientBasis{Polynomial{B,M}},FullBasis{B,M},SubBasis{B,M}}, + }, ) where {B,M} return Algebra{BT,B,M} end diff --git a/src/chebyshev.jl b/src/chebyshev.jl index eea1522..dc866a7 100644 --- a/src/chebyshev.jl +++ b/src/chebyshev.jl @@ -34,13 +34,15 @@ function (::Mul{Chebyshev})(a::MP.AbstractMonomial, b::MP.AbstractMonomial) vars_b = MP.variables(b) var_state_b = iterate(vars_b) while !isnothing(var_state_a) || !isnothing(var_state_b) - if isnothing(var_state_a) || (!isnothing(var_state_b) && var_state_b[1] > var_state_a[1]) + if isnothing(var_state_a) || + (!isnothing(var_state_b) && var_state_b[1] > var_state_a[1]) var_b, state_b = var_state_b for i in eachindex(terms) terms[i] = MA.mul!!(terms[i], var_b^MP.degree(b, var_b)) end var_state_b = iterate(vars_b, state_b) - elseif isnothing(var_state_b) || (!isnothing(var_state_a) && var_state_a[1] > var_state_b[1]) + elseif isnothing(var_state_b) || + (!isnothing(var_state_a) && var_state_a[1] > var_state_b[1]) var_a, state_a = var_state_a for i in eachindex(terms) terms[i] = MA.mul!!(terms[i], var_a^MP.degree(a, var_a)) @@ -67,7 +69,7 @@ function (::Mul{Chebyshev})(a::MP.AbstractMonomial, b::MP.AbstractMonomial) end function _add_mul_scalar_vector!(res, ::SubBasis, scalar, vector) - MA.operate!(MA.add_mul, res, scalar, vector) + return MA.operate!(MA.add_mul, res, scalar, vector) end function _add_mul_scalar_vector!(res, ::FullBasis, scalar, vector) @@ -76,7 +78,12 @@ function _add_mul_scalar_vector!(res, ::FullBasis, scalar, vector) end end -function SA.coeffs!(res, cfs, source::MonomialIndexedBasis{Chebyshev}, target::MonomialIndexedBasis{Monomial}) +function SA.coeffs!( + res, + cfs, + source::MonomialIndexedBasis{Chebyshev}, + target::MonomialIndexedBasis{Monomial}, +) MA.operate!(zero, res) for (k, v) in SA.nonzero_pairs(cfs) _add_mul_scalar_vector!(res, target, v, SA.coeffs(source[k], target)) @@ -85,7 +92,11 @@ function SA.coeffs!(res, cfs, source::MonomialIndexedBasis{Chebyshev}, target::M return res end -function SA.coeffs(cfs, source::SubBasis{Monomial}, target::FullBasis{Chebyshev}) +function SA.coeffs( + cfs, + source::SubBasis{Monomial}, + target::FullBasis{Chebyshev}, +) sub = explicit_basis_covering(target, source) # Need to make A square so that it's UpperTriangular extended = SubBasis{Monomial}(sub.monomials) @@ -94,11 +105,22 @@ function SA.coeffs(cfs, source::SubBasis{Monomial}, target::FullBasis{Chebyshev} 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) + return SA.SparseCoefficients( + sub.monomials, + LinearAlgebra.UpperTriangular(A) \ ext, + ) end -function SA.coeffs(cfs, source::FullBasis{Monomial}, target::FullBasis{Chebyshev}) - return SA.coeffs(SA.values(cfs), SubBasis{Monomial}(collect(SA.keys(cfs))), target) +function SA.coeffs( + cfs, + source::FullBasis{Monomial}, + target::FullBasis{Chebyshev}, +) + return SA.coeffs( + SA.values(cfs), + SubBasis{Monomial}(collect(SA.keys(cfs))), + target, + ) end function degree_one_univariate_polynomial( diff --git a/src/monomial.jl b/src/monomial.jl index be1f152..eee8e2c 100644 --- a/src/monomial.jl +++ b/src/monomial.jl @@ -125,7 +125,10 @@ end implicit_basis(::SubBasis{B,M}) where {B,M} = FullBasis{B,M}() implicit_basis(basis::FullBasis) = basis -function MA.promote_operation(::typeof(implicit_basis), ::Type{<:Union{FullBasis{B,M},SubBasis{B,M}}}) where {B,M} +function MA.promote_operation( + ::typeof(implicit_basis), + ::Type{<:Union{FullBasis{B,M},SubBasis{B,M}}}, +) where {B,M} return FullBasis{B,M} end @@ -135,7 +138,9 @@ end _explicit_basis(_, basis::SubBasis) = basis -explicit_basis(a::SA.AlgebraElement) = _explicit_basis(SA.coeffs(a), SA.basis(a)) +function explicit_basis(a::SA.AlgebraElement) + return _explicit_basis(SA.coeffs(a), SA.basis(a)) +end function explicit_basis_type(::Type{FullBasis{B,M}}) where {B,M} return SubBasis{B,M,MP.monomial_vector_type(M)} @@ -171,31 +176,43 @@ function MA.promote_operation( ) where {P<:MP.AbstractPolynomialLike} M = MP.monomial_type(P) T = MP.coefficient_type(P) - return SA.SparseCoefficients{ - M, - T, - MP.monomial_vector_type(M), - Vector{T}, - } + return SA.SparseCoefficients{M,T,MP.monomial_vector_type(M),Vector{T}} end function algebra_element(p::MP.AbstractPolynomialLike) - return algebra_element(sparse_coefficients(p), FullBasis{Monomial,MP.monomial_type(p)}()) + return algebra_element( + sparse_coefficients(p), + FullBasis{Monomial,MP.monomial_type(p)}(), + ) end function algebra_element(p::Polynomial{B,M}) where {B,M} - return algebra_element(sparse_coefficients(MP.term(1, p.monomial)), FullBasis{B,M}()) + return algebra_element( + sparse_coefficients(MP.term(1, p.monomial)), + FullBasis{B,M}(), + ) end function algebra_element(f::Function, basis::SubBasis) return algebra_element(map(f, eachindex(basis)), basis) end -function constant_algebra_element(::Type{FullBasis{B,M}}, ::Type{T}) where {B,M,T} - return algebra_element(sparse_coefficients(MP.polynomial(MP.term(one(T), MP.constant_monomial(M)))), FullBasis{B,M}()) +function constant_algebra_element( + ::Type{FullBasis{B,M}}, + ::Type{T}, +) where {B,M,T} + return algebra_element( + sparse_coefficients( + MP.polynomial(MP.term(one(T), MP.constant_monomial(M))), + ), + FullBasis{B,M}(), + ) end -function constant_algebra_element(::Type{<:SubBasis{B,M}}, ::Type{T}) where {B,M,T} +function constant_algebra_element( + ::Type{<:SubBasis{B,M}}, + ::Type{T}, +) where {B,M,T} return algebra_element([one(T)], SubBasis{B}([MP.constant_monomial(M)])) end @@ -240,7 +257,9 @@ function basis_covering_monomials( return SubBasis{B}(monos) end -Base.adjoint(p::Polynomial{B}) where {B<:AbstractMonomialIndexed} = Polynomial{B}(adjoint(p.monomial)) +function Base.adjoint(p::Polynomial{B}) where {B<:AbstractMonomialIndexed} + return Polynomial{B}(adjoint(p.monomial)) +end """ struct Monomial <: AbstractMonomialIndexed end @@ -298,8 +317,10 @@ end function _assert_constant(α) end -function _assert_constant(x::Union{Polynomial,SA.AlgebraElement,MP.AbstractPolynomialLike}) - error("Expected constant element, got type `$(typeof(x))`") +function _assert_constant( + x::Union{Polynomial,SA.AlgebraElement,MP.AbstractPolynomialLike}, +) + return error("Expected constant element, got type `$(typeof(x))`") end #function MA.operate!(::SA.UnsafeAddMul{<:Mul{Monomial}}, p::MP.AbstractPolynomial, args::Vararg{Any,N}) where {N} @@ -329,7 +350,12 @@ end # return a #end -function MA.operate!(::SA.UnsafeAddMul{typeof(*)}, a::SA.AlgebraElement{<:Algebra{<:MonomialIndexedBasis,Monomial}}, α, x::Polynomial{Monomial}) +function MA.operate!( + ::SA.UnsafeAddMul{typeof(*)}, + a::SA.AlgebraElement{<:Algebra{<:MonomialIndexedBasis,Monomial}}, + α, + x::Polynomial{Monomial}, +) _assert_constant(α) SA.unsafe_push!(a, x, α) return a diff --git a/src/orthogonal.jl b/src/orthogonal.jl index a636a7c..df7c42e 100644 --- a/src/orthogonal.jl +++ b/src/orthogonal.jl @@ -178,14 +178,22 @@ function MP.coefficients( end end -function SA.coeffs(p::Polynomial{B}, basis::SubBasis{Monomial}) where {B<:AbstractMultipleOrthogonal} +function SA.coeffs( + p::Polynomial{B}, + basis::SubBasis{Monomial}, +) where {B<:AbstractMultipleOrthogonal} full = implicit_basis(basis) return SA.coeffs(algebra_element(SA.coeffs(p, full), full), basis) end -function SA.coeffs(p::Polynomial{B}, ::FullBasis{Monomial}) where {B<:AbstractMultipleOrthogonal} - return sparse_coefficients(prod( - univariate_orthogonal_basis(B, var, deg)[deg+1] for - (var, deg) in MP.powers(p.monomial) - )) +function SA.coeffs( + p::Polynomial{B}, + ::FullBasis{Monomial}, +) where {B<:AbstractMultipleOrthogonal} + return sparse_coefficients( + prod( + univariate_orthogonal_basis(B, var, deg)[deg+1] for + (var, deg) in MP.powers(p.monomial) + ), + ) end diff --git a/src/polynomial.jl b/src/polynomial.jl index 755e9b6..5c99d8d 100644 --- a/src/polynomial.jl +++ b/src/polynomial.jl @@ -20,11 +20,7 @@ function Base.iterate(t::MP.Term, state) return nothing end end -function SA.unsafe_push!( - p::MP.AbstractPolynomial, - mono::MP.AbstractMonomial, - α, -) +function SA.unsafe_push!(p::MP.AbstractPolynomial, mono::MP.AbstractMonomial, α) return MA.operate!(MA.add_mul, p, α, mono) end function MA.operate!( @@ -154,13 +150,21 @@ function MA.operate_to!( return p end -MP.polynomial(a::SA.AbstractCoefficients) = MP.polynomial(SA.values(a), SA.keys(a)) +function MP.polynomial(a::SA.AbstractCoefficients) + return MP.polynomial(SA.values(a), SA.keys(a)) +end function MP.polynomial(a::SA.AlgebraElement) - return MP.polynomial(SA.coeffs(a, FullBasis{Monomial,MP.monomial_type(typeof(a))}())) + return MP.polynomial( + SA.coeffs(a, FullBasis{Monomial,MP.monomial_type(typeof(a))}()), + ) end -function Base.isapprox(p::MP.AbstractPolynomialLike, a::SA.AlgebraElement; kws...) +function Base.isapprox( + p::MP.AbstractPolynomialLike, + a::SA.AlgebraElement; + kws..., +) return isapprox(p, MP.polynomial(a); kws...) end @@ -168,10 +172,18 @@ function Base.isapprox(a::SA.AlgebraElement, b::SA.AlgebraElement; kws...) return isapprox(MP.polynomial(a), b; kws...) end -function Base.isapprox(a::SA.AlgebraElement, p::MP.AbstractPolynomialLike; kws...) +function Base.isapprox( + a::SA.AlgebraElement, + p::MP.AbstractPolynomialLike; + kws..., +) return isapprox(p, a; kws...) end function Base.isapprox(a::SA.AlgebraElement, α::Number; kws...) - return isapprox(a, α * constant_algebra_element(typeof(SA.basis(a)), typeof(α)); kws...) + return isapprox( + a, + α * constant_algebra_element(typeof(SA.basis(a)), typeof(α)); + kws..., + ) end diff --git a/src/quotient.jl b/src/quotient.jl index 92b882d..2ac4c33 100644 --- a/src/quotient.jl +++ b/src/quotient.jl @@ -15,28 +15,51 @@ function MP.coefficients(p, basis::QuotientBasis) return MP.coefficients(rem(p, basis.divisor), basis.basis) end -function SA.coeffs(p, ::FullBasis{Monomial}, basis::Union{QuotientBasis,SubBasis}) +function SA.coeffs( + p, + ::FullBasis{Monomial}, + basis::Union{QuotientBasis,SubBasis}, +) return MP.coefficients(MP.polynomial(p), basis) end -function SA.coeffs(coeffs, sub::SubBasis{Monomial}, basis::Union{MonomialIndexedBasis,QuotientBasis}) +function SA.coeffs( + coeffs, + sub::SubBasis{Monomial}, + basis::Union{MonomialIndexedBasis,QuotientBasis}, +) return SA.coeffs(MP.polynomial(coeffs, sub.monomials), parent(sub), basis) end -function MA.operate!(op::SA.UnsafeAddMul{typeof(*)}, a::SA.AlgebraElement{<:Algebra{<:QuotientBasis,Monomial}}, α, p::Polynomial{Monomial}) +function MA.operate!( + op::SA.UnsafeAddMul{typeof(*)}, + a::SA.AlgebraElement{<:Algebra{<:QuotientBasis,Monomial}}, + α, + p::Polynomial{Monomial}, +) _assert_constant(α) for t in MP.terms(rem(p.monomial, SA.basis(a).divisor)) - MA.operate!(op, algebra_element(SA.coeffs(a), SA.basis(a).basis), α * MP.coefficient(t), Polynomial{Monomial}(MP.monomial(t))) + MA.operate!( + op, + algebra_element(SA.coeffs(a), SA.basis(a).basis), + α * MP.coefficient(t), + Polynomial{Monomial}(MP.monomial(t)), + ) #SA.unsafe_push!(SA.coeffs(a), SA.basis(a).basis[Polynomial{Monomial}(MP.monomial(k))], α * MP.coefficient(t)) end return a end -function SA.adjoint_coeffs(coeffs, src::SubBasis{Monomial}, dest::QuotientBasis{<:Polynomial{Monomial}}) +function SA.adjoint_coeffs( + coeffs, + src::SubBasis{Monomial}, + dest::QuotientBasis{<:Polynomial{Monomial}}, +) return map(src.monomials) do mono return sum( - MP.coefficient(t) * coeffs[dest.basis[Polynomial{Monomial}(MP.monomial(t))]] - for t in MP.terms(rem(mono, dest.divisor)) + MP.coefficient(t) * + coeffs[dest.basis[Polynomial{Monomial}(MP.monomial(t))]] for + t in MP.terms(rem(mono, dest.divisor)) ) end end diff --git a/src/scaled.jl b/src/scaled.jl index 4d5a95f..4e7e51f 100644 --- a/src/scaled.jl +++ b/src/scaled.jl @@ -35,8 +35,11 @@ struct ScaledMonomial <: AbstractMonomial end function (::Mul{ScaledMonomial})(a::MP.AbstractMonomial, b::MP.AbstractMonomial) mono = a * b - α = prod(MP.variables(mono); init = inv(binomial(MP.degree(mono), MP.degree(a)))) do v - binomial(MP.degree(mono, v), MP.degree(a, v)) + α = prod( + MP.variables(mono); + init = inv(binomial(MP.degree(mono), MP.degree(a))), + ) do v + return binomial(MP.degree(mono, v), MP.degree(a, v)) end return sparse_coefficients(MP.term(√α, mono)) end @@ -69,7 +72,11 @@ function scaling(m::MP.AbstractMonomial) return √(factorial(MP.degree(m)) / prod(factorial, MP.exponents(m))) end unscale_coef(t::MP.AbstractTerm) = MP.coefficient(t) / scaling(MP.monomial(t)) -function SA.coeffs(t::MP.AbstractTermLike, ::FullBasis{ScaledMonomial}, ::FullBasis{Monomial}) +function SA.coeffs( + t::MP.AbstractTermLike, + ::FullBasis{ScaledMonomial}, + ::FullBasis{Monomial}, +) mono = MP.monomial(t) return MP.term(mono * MP.coefficient(t), mono) end @@ -79,25 +86,47 @@ end function MP.coefficients(p, basis::SubBasis{ScaledMonomial}) return MP.coefficients(p, basis.monomials) ./ scaling.(MP.monomials(p)) end -function SA.coeffs(p::MP.AbstractPolynomialLike, ::FullBasis{Monomial}, basis::FullBasis{ScaledMonomial}) +function SA.coeffs( + p::MP.AbstractPolynomialLike, + ::FullBasis{Monomial}, + basis::FullBasis{ScaledMonomial}, +) return MP.polynomial(MP.coefficients(p, basis), MP.monomials(p)) end -function SA.coeffs!(res, cfs, source::MonomialIndexedBasis{ScaledMonomial}, target::MonomialIndexedBasis{Monomial}) +function SA.coeffs!( + res, + cfs, + source::MonomialIndexedBasis{ScaledMonomial}, + target::MonomialIndexedBasis{Monomial}, +) MA.operate!(zero, res) for (k, v) in SA.nonzero_pairs(cfs) mono = source[k].monomial - SA.unsafe_push!(res, target[Polynomial{Monomial}(mono)], v * scaling(mono)) + SA.unsafe_push!( + res, + target[Polynomial{Monomial}(mono)], + v * scaling(mono), + ) end MA.operate!(SA.canonical, res) return res end -function SA.coeffs!(res, cfs, source::MonomialIndexedBasis{Monomial}, target::MonomialIndexedBasis{ScaledMonomial}) +function SA.coeffs!( + res, + cfs, + source::MonomialIndexedBasis{Monomial}, + target::MonomialIndexedBasis{ScaledMonomial}, +) MA.operate!(zero, res) for (k, v) in SA.nonzero_pairs(cfs) mono = source[k].monomial - SA.unsafe_push!(res, target[Polynomial{ScaledMonomial}(mono)], v / scaling(mono)) + SA.unsafe_push!( + res, + target[Polynomial{ScaledMonomial}(mono)], + v / scaling(mono), + ) end MA.operate!(SA.canonical, res) return res diff --git a/test/chebyshev.jl b/test/chebyshev.jl index c66542f..a560624 100644 --- a/test/chebyshev.jl +++ b/test/chebyshev.jl @@ -9,7 +9,8 @@ using DynamicPolynomials @test b.coeffs == 1 // 2 + 1 // 2 * x^2 c = b * b @test c.coeffs == 3 // 8 + 1 // 2 * x^2 + 1 // 8 * x^4 - @test a * MB.Polynomial{MB.Chebyshev}(constant_monomial(typeof(x))) == a * MB.Polynomial{MB.Chebyshev}(x^0) + @test a * MB.Polynomial{MB.Chebyshev}(constant_monomial(typeof(x))) == + a * MB.Polynomial{MB.Chebyshev}(x^0) end @testset "Orthogonal" begin