Skip to content

Commit

Permalink
Fix format
Browse files Browse the repository at this point in the history
  • Loading branch information
blegat committed Jun 14, 2024
1 parent a105c3e commit fee4ab8
Show file tree
Hide file tree
Showing 8 changed files with 186 additions and 64 deletions.
15 changes: 8 additions & 7 deletions src/MultivariateBases.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
38 changes: 30 additions & 8 deletions src/chebyshev.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand All @@ -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)
Expand All @@ -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))
Expand All @@ -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)
Expand All @@ -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(
Expand Down
60 changes: 43 additions & 17 deletions src/monomial.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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)}
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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}
Expand Down Expand Up @@ -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
Expand Down
20 changes: 14 additions & 6 deletions src/orthogonal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
32 changes: 22 additions & 10 deletions src/polynomial.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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!(
Expand Down Expand Up @@ -154,24 +150,40 @@ 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

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
37 changes: 30 additions & 7 deletions src/quotient.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Loading

0 comments on commit fee4ab8

Please sign in to comment.