improve computing with expansions and basis functions

Project.toml

@@ -1,7 +1,7 @@
 name = "BasisFunctions"
 uuid = "4343a256-5453-507d-8aad-01a9d7189916"
 authors = ["Daan Huybrechs <daan.huybrechs@kuleuven.be>"]
-version = "0.6.6"
+version = "0.7"
 
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
@@ -40,7 +40,7 @@ BasisFunctionsPGFPlotsXExt = "PGFPlotsX"
 BlockArrays = "0.16"
 Calculus = "0.5"
 CompositeTypes = "0.1.3"
-DomainIntegrals = "0.4.5"
+DomainIntegrals = "0.5"
 DomainSets = "0.7"
 DoubleFloats = "1"
 FFTW = "1.6"
@@ -63,7 +63,7 @@ SpecialFunctions = "1.0, 2.0"
 StaticArrays = "1.4"
 Test = "1"
 ToeplitzMatrices = "0.7"
-julia = "1.9"
+julia = "1.10"
 
 [extras]
 Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"

src/BasisFunctions.jl

@@ -56,7 +56,6 @@ import LinearAlgebra:
     norm, pinv, normalize, cross, ×, dot, adjoint, mul!, rank,
     diag, isdiag, eigvals, issymmetric, svdvals
 import SparseArrays: sparse
-@deprecate sparse_matrix sparse
 
 
 import BlockArrays.BlockVector
@@ -158,6 +157,7 @@ export GenericFunctionSpace, space, MeasureSpace, FourierSpace, ChebyshevTSpace,
 # from bases/dictionary/dictionary.jl
 export Dictionary, Dictionary1d, Dictionary2d, Dictionary3d,
     interpolation_grid, left, right, domain, codomain,
+    domaintype, codomaintype, coefficienttype,
     measure, hasmeasure,
     eval_expansion, eval_element, eval_element_derivative, eval_gradient,
     resize,
@@ -356,8 +356,13 @@ include("bases/dictionary/expansions.jl")
 
 include("util/products.jl")
 
-include("computations/generic_operators.jl")
+include("computations/transform.jl")
+include("computations/evaluation.jl")
+include("computations/approximation.jl")
+include("computations/differentiation.jl")
 include("computations/conversion.jl")
+include("computations/generic_operators.jl")
+include("computations/fun.jl")
 
 ################################################################
 # Generic dictionaries
@@ -418,4 +423,7 @@ include("util/plot.jl")
 include("test/Test.jl")
 
 include("examples/pwconstants.jl")
+
+include("deprecated.jl")
+
 end # module

src/bases/dictionary/dict_evaluation.jl

@@ -217,8 +217,6 @@ default_eval_expansion(dict::Dictionary, coefficients, x) =
 default_unsafe_eval_expansion(dict::Dictionary, coefficients, x) =
     sum(coefficients[idx]*val for (idx,val) in unsafe_pointvalues(dict, x))
 
-@deprecate eval_expansion(dict::Dictionary, coefficients, grid::AbstractGrid; options...) eval_expansion.(Ref(dict), Ref(coefficients), grid; options...)
-
 function eval_expansion_derivative(dict::Dictionary, coefficients, x, order; options...)
     @assert size(coefficients) == size(dict)
     in_support(dict, x) ? unsafe_eval_expansion_derivative(dict, coefficients, x, order) : zero(span_codomaintype(dict, coefficients))

src/bases/dictionary/dictionary.jl

@@ -1,10 +1,4 @@
 
-export domaintype,
-    codomaintype,
-    coefficienttype,
-    prectype
-
-
 """
 A `Dictionary{S,T}` is an ordered family of functions, in which each function
 maps a variable of type `S` to a variable of type `T`. The dictionary can be
@@ -101,6 +95,11 @@ lastindex(d::Dictionary) = last(eachindex(d))
 
 ≈(d1::Dictionary, d2::Dictionary) = d1==d2
 
+isequaldict(d1, d2) = isequaldict1(d1, d2)
+isequaldict1(d1, d2) = isequaldict2(d1, d2)
+isequaldict2(d1, d2) = default_isequaldict(d1, d2)
+default_isequaldict(d1, d2) = d1 == d2
+
 #############################
 # Domain and codomain type
 #############################

src/bases/dictionary/expansions.jl

@@ -54,6 +54,7 @@ function expansion(dict::Dictionary, coefficients)
     T = promote_type(coefficienttype(dict), eltype(coef))
     Expansion(ensure_coefficienttype(T, dict), coef)
 end
+expansion(e::Expansion) = e
 
 similar(e::Expansion, coefficients) = Expansion(dictionary(e), coefficients)
 
@@ -64,7 +65,6 @@ coefficients(e::Expansion) = e.coefficients
 mapped_expansion(e::Expansion, m) =
     Expansion(mapped_dict(dictionary(e), m), coefficients(e))
 
-export random_expansion
 random_expansion(dict::Dictionary) = Expansion(dict, rand(dict))
 
 # For expansions of composite types, return a Expansion of a subdict
@@ -119,17 +119,6 @@ ei(dim,i, coefficients) = tuple((coefficients*Matrix{Int}(I, dim, dim)[:,i])...)
 differentiate(f::Expansion, var, order) = differentiate(f, ei(dimension(f), var, order))
 antidifferentiate(f::Expansion, var, order) = antidifferentiate(f, ei(dimension(f), var, order))
 
-function Base.:^(f::Expansion, i::Int)
-    @assert i >= 0
-    if i == 1
-        f
-    elseif i == 2
-        f * f
-    else
-        f^(i-1) * f
-    end
-end
-
 # To be implemented: Laplacian (needs multiplying functions)
 ## Δ(f::Expansion)
 
@@ -140,8 +129,6 @@ adjoint(f::Expansion) = differentiate(f)
 
 roots(f::Expansion) = expansion_roots(dictionary(f), coefficients(f))
 
-Base.@deprecate roots(dict::Dictionary, coef::AbstractVector) expansion_roots(dict, coef)
-
 show(io::IO, mime::MIME"text/plain", fun::Expansion) = composite_show(io, mime, fun)
 Display.displaystencil(fun::Expansion) = ["Expansion(", dictionary(fun), ", ", coefficients(fun), ")"]
 
@@ -199,56 +186,12 @@ function to_compatible_expansions(e1::Expansion, e2::Expansion)
     end
 end
 
-function (+)(f::Expansion, g::Expansion)
-    dict, coef = expansion_sum(dictionary(f),dictionary(g),coefficients(f),coefficients(g))
-    Expansion(dict, coef)
-end
-(-)(e1::Expansion, e2::Expansion) = e1 + similar(e2, -coefficients(e2))
-
-(-)(e::Expansion) = similar(e, -coefficients(e))
-
-@deprecate dict_multiply expansion_multiply
-expansion_multiply(dict1, dict2, coef1, coef2) =
-    expansion_multiply1(dict1, dict2, coef1, coef2)
-expansion_multiply1(dict1, dict2, coef1, coef2) =
-    expansion_multiply2(dict1, dict2, coef1, coef2)
-expansion_multiply2(dict1, dict2, coef1, coef2) =
-    error("Multiplication of expansions in these dictionaries is not implemented.")
-
-expansion_sum(dict1, dict2, coef1, coef2) =
-    expansion_sum1(dict1, dict2, coef1, coef2)
-expansion_sum1(dict1, dict2, coef1, coef2) =
-    expansion_sum2(dict1, dict2, coef1, coef2)
-function expansion_sum2(dict1, dict2, coef1, coef2)
-    if iscompatible(dict1, dict2)
-        default_expansion_sum(dict1, dict2, coef1, coef2)
-    else
-        error("Summation of expansions in these dictionaries is not implemented.")
-    end
-end
-
-function default_expansion_sum(dict1, dict2, coef1, coef2)
-    @assert iscompatible(dict1, dict2)
-    e1, e2 = to_compatible_expansions(Expansion(dict1, coef1), Expansion(dict2, coef2))
-    dictionary(e1), coefficients(e1)+coefficients(e2)
-end
-
-function (*)(f::Expansion, g::Expansion)
-    dict, coef = expansion_multiply(dictionary(f),dictionary(g),coefficients(f),coefficients(g))
-    Expansion(dict, coef)
-end
-
 (*)(op::DictionaryOperator, e::Expansion) = apply(op, e)
 
 Base.complex(e::Expansion) = Expansion(complex(dictionary(e)), complex(coefficients(e)))
-function Base.real(e::Expansion)
-    dict, coef = expansion_real(dictionary(e), coefficients(e))
-    Expansion(dict, coef)
-end
-function Base.imag(e::Expansion)
-    dict, coef = expansion_imag(dictionary(e), coefficients(e))
-    Expansion(dict, coef)
-end
+Base.real(e::Expansion) = Expansion(expansion_real(dictionary(e), coefficients(e))...)
+Base.imag(e::Expansion) = Expansion(expansion_imag(dictionary(e), coefficients(e))...)
+
 function expansion_real(dict::Dictionary, coefficients)
     @assert isreal(dict)
     real(dict), real(coefficients)
@@ -263,14 +206,6 @@ Base.one(e::Expansion) = similar(e, coefficients_of_one(dictionary(e)))
 expansion_of_one(dict::Dictionary) = Expansion(dict, coefficients_of_one(dict))
 expansion_of_x(dict::Dictionary) = Expansion(dict, coefficients_of_x(dict))
 
-for op in (:+, :-)
-    @eval Base.$op(a::Number, e::Expansion) = $op(a*one(e), e)
-    @eval Base.$op(e::Expansion, a::Number) = $op(e, a*one(e))
-end
-Base.:*(e::Expansion, a::Number) = expansion(dictionary(e), coefficients(e)*a)
-Base.:*(a::Number, e::Expansion) = expansion(dictionary(e), a*coefficients(e))
-Base.:/(e::Expansion, a::Number) = expansion(dictionary(e), coefficients(e)/a)
-
 apply(op::DictionaryOperator, e::Expansion) = Expansion(dest(op), op * coefficients(e))
 
 # TODO: we may not want to identify an Expansion with its array of coefficients
@@ -281,22 +216,3 @@ Base.collect(e::Expansion) = coefficients(e)
 
 Base.BroadcastStyle(e::Expansion) = Base.Broadcast.DefaultArrayStyle{dimension(e)}()
 
-# Interoperate with basis functions.
-# TODO: implement cleaner using promotion mechanism
-
-Base.:*(φ1::AbstractBasisFunction, φ2::AbstractBasisFunction) = expansion(φ1) * expansion(φ2)
-Base.:*(f::Expansion, φ2::AbstractBasisFunction) = f * expansion(φ2)
-Base.:*(φ1::AbstractBasisFunction, f::Expansion) = expansion(φ1) * f
-
-Base.:*(A::Number, φ::AbstractBasisFunction) = A * expansion(φ)
-Base.:*(φ::AbstractBasisFunction, A::Number) = A * expansion(φ)
-
-Base.:+(φ1::AbstractBasisFunction, φ2::AbstractBasisFunction) = expansion(φ1) + expansion(φ2)
-Base.:+(f::Expansion, φ2::AbstractBasisFunction) = f + expansion(φ2)
-Base.:+(φ1::AbstractBasisFunction, f::Expansion) = expansion(φ1) + f
-
-Base.:-(φ1::AbstractBasisFunction, φ2::AbstractBasisFunction) = expansion(φ1) - expansion(φ2)
-Base.:-(f::Expansion, φ2::AbstractBasisFunction) = f - expansion(φ2)
-Base.:-(φ1::AbstractBasisFunction, f::Expansion) = expansion(φ1) - f
-
-

src/bases/dictionary/gram.jl

@@ -8,12 +8,6 @@ hasmeasure(Φ::Dictionary) = false
 # Shortcut: Dictionaries of the same type have just one measure
 defaultmeasure(Φ1::D, Φ2::D) where {D <: Dictionary} = measure(Φ1)
 
-@deprecate innerproduct(Φ1::Dictionary, i, Φ2::Dictionary, j; options...) dict_innerproduct(Φ1, i, Φ2, j; options...)
-@deprecate innerproduct(Φ1::Dictionary, i, Φ2::Dictionary, j, measure; options...) dict_innerproduct(Φ1, i, Φ2, j, measure; options...)
-@deprecate innerproduct_native(Φ1::Dictionary, i, Φ2::Dictionary, j, measure; options...) dict_innerproduct_native(Φ1, i, Φ2, j, measure; options...)
-@deprecate innerproduct1(Φ1::Dictionary, i, Φ2::Dictionary, j, measure; options...) dict_innerproduct1(Φ1, i, Φ2, j, measure; options...)
-@deprecate innerproduct2(Φ1::Dictionary, i, Φ2::Dictionary, j, measure; options...) dict_innerproduct2(Φ1, i, Φ2, j, measure; options...)
-
 dict_innerproduct(Φ1::Dictionary, i, Φ2::Dictionary, j; options...) =
     dict_innerproduct(Φ1, i, Φ2, j, defaultmeasure(Φ1, Φ2); options...)
 

src/bases/dictionary/promotion.jl

@@ -90,3 +90,20 @@ widen_prectype(::Type{Tuple{A,B}}, ::Type{U}) where {A,B,U} =
     Tuple{widen_prectype(A,U),widen_prectype(B,U)}
 widen_prectype(::Type{Tuple{A,B,C}}, ::Type{U}) where {A,B,C,U} =
     Tuple{widen_prectype(A,U),widen_prectype(B,U),widen_prectype(C,U)}
+
+
+## Promotion with conversion
+
+"""
+Promote two dictionaries to a common family, in such a way that expansions in both given
+dictionaries can be converted to expansions in the promoted dictionaries.
+"""
+promote_convertible(d1, d2) = promote_convertible1(promote_domaintype(d1, d2)...)
+promote_convertible1(d1, d2) = promote_convertible2(d1, d2)
+promote_convertible2(d1, d2) = default_promote_convertible(d1, d2)
+
+nopromotion(d1,d2) = error("Don't know how to promote and convert between dictionaries $(d1) and $(d2).")
+
+default_promote_convertible(d1::D, d2::D) where D = d1, d2
+default_promote_convertible(d1, d2) =
+    isequaldict(d1,d2) ? d1 : nopromotion(d1,d2)

src/bases/fourier/cosineseries.jl

@@ -17,6 +17,7 @@ show(io::IO, b::CosineSeries{Float64}) = print(io, "CosineSeries($(length(b)))")
 show(io::IO, b::CosineSeries{T}) where T = print(io, "CosineSeries{$(T)}($(length(b)))")
 
 similar(b::CosineSeries, ::Type{T}, n::Int) where {T} = CosineSeries{T}(n)
+isequaldict(b1::CosineSeries, b2::CosineSeries) = length(b1)==length(b2)
 
 isbasis(b::CosineSeries) = true
 isorthogonal(b::CosineSeries, ::FourierWeight) = true

src/bases/fourier/fourier.jl

@@ -93,6 +93,7 @@ tofourier(d::FourierLike{T}) where {T} = convert(Fourier{T}, d)
 size(b::Fourier) = (b.n,)
 
 similar(b::Fourier, ::Type{T}, n::Int) where {T} = Fourier{T}(n)
+isequaldict(b1::Fourier, b2::Fourier) = length(b1)==length(b2)
 
 show(io::IO, b::Fourier{Float64}) = print(io, "Fourier($(length(b)))")
 show(io::IO, b::Fourier{T}) where T = print(io, "Fourier{$(T)}($(length(b)))")

src/bases/fourier/sinc.jl

@@ -16,6 +16,7 @@ show(io::IO, b::PeriodicSincFunctions{T}) where T =
     print(io, "PeriodicSincFunctions{$(T)}($(length(b)))")
 
 similar(b::PeriodicSincFunctions, ::Type{T}, n::Int) where {T} = PeriodicSincFunctions{T}(n)
+isequaldict(b1::PeriodicSincFunctions, b2::PeriodicSincFunctions) = length(b1)==length(b2)
 
 size(b::PeriodicSincFunctions) = (b.n,)
 

src/bases/fourier/sineseries.jl

@@ -14,6 +14,7 @@ show(io::IO, b::SineSeries{T}) where T = print(io, "SineSeries{$(T)}($(length(b)
 SineSeries(n::Int) = SineSeries{Float64}(n)
 
 similar(b::SineSeries, ::Type{T}, n::Int) where {T} = SineSeries{T}(n)
+isequaldict(b1::SineSeries, b2::SineSeries) = length(b1)==length(b2)
 
 isbasis(b::SineSeries) = true
 isorthogonal(b::SineSeries, ::FourierWeight) = true

src/bases/fourier/trig.jl

@@ -21,6 +21,7 @@ coshalf(n::Int) = 1+nhalf(n)
 sinhalf(n::Int) = n-coshalf(n)
 
 similar(d::TrigSeries, ::Type{T}, n::Int) where {T} = TrigSeries{T}(n)
+isequaldict(b1::TrigSeries, b2::TrigSeries) = length(b1)==length(b2)
 
 show(io::IO, d::TrigSeries{Float64}) = print(io, "TrigSeries($(length(d)))")
 show(io::IO, d::TrigSeries{T}) where T = print(io, "TrigSeries{$(T)}($(length(d)))")

src/bases/modifiers/complexified_dict.jl

@@ -22,6 +22,9 @@ Base.real(dict::ComplexifiedDict) = superdict(dict)
 
 similardictionary(s::ComplexifiedDict, s2::Dictionary) = ComplexifiedDict(s2)
 
+isequaldict1(d1::ComplexifiedDict, d2) = isequaldict(superdict(d1), d2)
+isequaldict2(d1, d2::ComplexifiedDict) = isequaldict(d1, superdict(d2))
+
 components(dict::ComplexifiedDict) = map(ComplexifiedDict, components(superdict(dict)))
 component(dict::ComplexifiedDict, i) = ComplexifiedDict(component(superdict(dict), i))
 
@@ -58,6 +61,25 @@ span_issubset2(d1, d2::ComplexifiedDict) = span_issubset(d1, superdict(d2))
 iscompatible(d1::Dictionary, d2::ComplexifiedDict) = iscompatible(d1, superdict(d2))
 iscompatible(d1::ComplexifiedDict, d2::Dictionary) = iscompatible(superdict(d1), d2)
 iscompatible(d1::ComplexifiedDict, d2::ComplexifiedDict) = iscompatible(superdict(d1), superdict(d2))
+
+function promote_convertible1(d1::ComplexifiedDict, d2::ComplexifiedDict)
+    reald1, reald2 = promote_convertible(superdict(d1), superdict(d2))
+    complex(reald1), complex(reald2)
+end
+function promote_convertible1(d1::ComplexifiedDict, d2)
+    reald1, reald2 = promote_convertible(superdict(d1), d2)
+    complex(reald1), complex(reald2)
+end
+function promote_convertible2(d1, d2::ComplexifiedDict)
+    reald1, reald2 = promote_convertible(d1, superdict(d2))
+    complex(reald1), complex(reald2)
+end
+
+function conversion1(T, d1::ComplexifiedDict, d2; options...)
+    @assert !isreal(T)
+    op = conversion(T, superdict(d1), d2; options...)
+    wrap_operator(d1, d2, op)
+end
 function conversion2(T, d1, d2::ComplexifiedDict; options...)
     @assert !isreal(T)
     op = conversion(T, d1, superdict(d2); options...)

src/bases/modifiers/mapped_dict.jl

@@ -25,9 +25,6 @@ mapped_dict2(dict, map) = MappedDict(dict, map)
 mapped_dict1(dict::MappedDict, map) = mapped_dict(superdict(dict), map ∘ forward_map(dict))
 mapped_dict2(dict, map::IdentityMap) = dict
 
-# Convenience function, similar to apply_map for grids etcetera
-@deprecate apply_map(dict::Dictionary, map) mapped_dict(dict, map)
-
 components(dict::MappedDict) = map(d->mapped_dict(d, forward_map(dict)), components(superdict(dict)))
 
 # In the constructor we check the domain and codomain types.

src/bases/modifiers/operated_dict.jl

@@ -283,6 +283,9 @@ function mixedgram2(T, dict1, dict2::OperatedDict, measure; options...)
 	wrap_operator(dict2, dict1, G *  operator(dict2))
 end
 
+promote_convertible1(d1::OperatedDict, d2) = promote_convertible(superdict(d1), d2)
+promote_convertible2(d1, d2::OperatedDict) = promote_convertible(d1, superdict(d2))
+
 function conversion1(::Type{T}, src::OperatedDict, dest; options...) where {T}
 	op = conversion(T, superdict(src), dest; options...) * operator(src)
 	wrap_operator(src, dest, op)

src/bases/modifiers/subdicts.jl

@@ -99,6 +99,20 @@ dict_innerproduct1(d1::Subdictionary, i, d2, j, measure; options...) =
 dict_innerproduct2(d1, i, d2::Subdictionary, j, measure; options...) =
     dict_innerproduct(d1, i, superdict(d2), superindices(d2, j), measure; options...)
 
+promote_convertible1(d1::Subdictionary, d2) =
+    d1 == d2 ? (d1, d2) : promote_convertible(superdict(d1), d2)
+promote_convertible2(d1, d2::Subdictionary) =
+    d1 == d2 ? (d1, d2) : promote_convertible(d1, superdict(d2))
+
+function conversion1(::Type{T}, d1::Subdictionary, d2) where T
+    if d1 == d2
+        IdentityOperator{T}(d1, d1)
+    else
+        E = extension(T, d1, superdict(d1))
+        op = conversion(superdict(d1), d2)
+        op * E
+    end
+end
 
 
 #####################