Derivative -> diff (#144)

* Derivative -> diff

* Derivative -> diff

* mapped Weighted diff

* Use grammatrix

* normalized diff

* Update jacobi.jl

* Lanczos mul

* Update interlace.jl

* Mapped weightedgrammatrix

* Add test from SingularIntegrals.jl

* fixes

* Update ultraspherical.jl

* Update ClassicalOrthogonalPolynomials.jl

* increase coverage

* Update test_legendre.jl

* Update legendre.jl

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.10.1"
+version = "0.11"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -29,7 +29,7 @@ ArrayLayouts = "1.0.1"
 BandedMatrices = "0.17.17"
 BlockArrays = "0.16.9"
 BlockBandedMatrices = "0.12"
-ContinuumArrays = "0.13"
+ContinuumArrays = "0.14"
 DomainSets = "0.6"
 FFTW = "1.1"
 FastGaussQuadrature = "0.5"
@@ -41,7 +41,7 @@ InfiniteLinearAlgebra = "0.6.16"
 IntervalSets = "0.7"
 LazyArrays = "1.0.1"
 LazyBandedMatrices = "0.8.5"
-QuasiArrays = "0.10"
+QuasiArrays = "0.11"
 SpecialFunctions = "1.0, 2"
 julia = "1.9"
 

examples/idealfluidflow.jl

@@ -11,7 +11,7 @@ import Base: @_inline_meta, axes, getindex, unsafe_getindex, convert, prod, *, /
                 IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy, setindex,
                 first, last, Slice, size, length, axes, IdentityUnitRange, sum, _sum, cumsum,
                 to_indices, tail, getproperty, inv, show, isapprox, summary,
-                findall, searchsortedfirst
+                findall, searchsortedfirst, diff
 import Base.Broadcast: materialize, BroadcastStyle, broadcasted, Broadcasted
 import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, adjointlayout,
                 sub_materialize, arguments, sub_paddeddata, paddeddata, PaddedLayout, resizedata!, LazyVector, ApplyLayout, call,
@@ -30,15 +30,15 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
                     LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle,
                     _getindex, layout_getindex, _factorize, AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector,
-                    AbstractQuasiFill, _dot, _equals, QuasiArrayLayout, PolynomialLayout
+                    AbstractQuasiFill, _equals, QuasiArrayLayout, PolynomialLayout, diff_layout
 
 import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, InfiniteCardinal, InfRanges
 import InfiniteLinearAlgebra: chop!, chop, pad, choplength, compatible_resize!
 import ContinuumArrays: Basis, Weight, basis_axes, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization,
     grid, plotgrid, _plotgrid, _grid, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, ExpansionLayout, basismap,
-    AffineQuasiVector, AffineMap, WeightLayout, AbstractWeightedBasisLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
+    AffineQuasiVector, AffineMap, AbstractWeightLayout, AbstractWeightedBasisLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
     checkpoints, weight, unweighted, MappedBasisLayouts, sum_layout, invmap, plan_ldiv, layout_broadcasted, MappedBasisLayout, SubBasisLayout, _broadcastbasis,
-    plan_transform, plan_grid_transform, MAX_PLOT_POINTS, MulPlan
+    plan_transform, plan_grid_transform, MAX_PLOT_POINTS, MulPlan, grammatrix, AdjointBasisLayout, grammatrix_layout
 import FastTransforms: Λ, forwardrecurrence, forwardrecurrence!, _forwardrecurrence!, clenshaw, clenshaw!,
                         _forwardrecurrence_next, _clenshaw_next, check_clenshaw_recurrences, ChebyshevGrid, chebyshevpoints, Plan, ScaledPlan, th_cheb2leg
 
@@ -160,7 +160,7 @@ end
 gives the singularity structure of an expansion, e.g.,
 `JacobiWeight`.
 """
-singularities(::WeightLayout, w) = w
+singularities(::AbstractWeightLayout, w) = w
 singularities(lay::BroadcastLayout, a::AbstractQuasiVector) = singularitiesbroadcast(call(a), map(singularities, arguments(lay, a))...)
 singularities(::WeightedBasisLayouts, a) = singularities(BroadcastLayout{typeof(*)}(), a)
 singularities(::WeightedOPLayout, a) = singularities(weight(a))
@@ -182,17 +182,22 @@ singularities(r::Base.RefValue) = r[] # pass through
 orthogonalityweight(P::SubQuasiArray{<:Any,2,<:Any,<:Tuple{AbstractAffineQuasiVector,Slice}}) =
     orthogonalityweight(parent(P))[parentindices(P)[1]]
 
-function massmatrix(P::SubQuasiArray{<:Any,2,<:Any,<:Tuple{AbstractAffineQuasiVector,Slice}})
+
+weighted(P::AbstractQuasiMatrix) = Weighted(P)
+
+weightedgrammatrix(P) = weightedgrammatrix_layout(MemoryLayout(P), P)
+function weightedgrammatrix_layout(::MappedOPLayout, P)
     Q = parent(P)
     kr,jr = parentindices(P)
-    massmatrix(Q)/kr.A
+    @assert kr isa AbstractAffineQuasiVector
+    weightedgrammatrix(Q)/kr.A
 end
 
-weighted(P::AbstractQuasiMatrix) = Weighted(P)
+grammatrix_layout(::MappedOPLayout, P) = grammatrix_layout(MappedBasisLayout(), P)
 
 OrthogonalPolynomial(w::Weight) =error("Override for $(typeof(w))")
 
-@simplify *(B::Identity, C::OrthogonalPolynomial) = C*jacobimatrix(C)
+@simplify *(B::Identity, C::OrthogonalPolynomial) = ApplyQuasiMatrix(*, C, jacobimatrix(C))
 
 function layout_broadcasted(::Tuple{PolynomialLayout,AbstractOPLayout}, ::typeof(*), x::Inclusion, C)
     x == axes(C,1) || throw(DimensionMismatch())

src/ClassicalOrthogonalPolynomials.jl

@@ -9,6 +9,10 @@ struct ChebyshevWeight{kind,T} <: AbstractJacobiWeight{T} end
 ChebyshevWeight{kind}() where kind = ChebyshevWeight{kind,Float64}()
 ChebyshevWeight() = ChebyshevWeight{1,Float64}()
 
+AbstractQuasiArray{T}(::ChebyshevWeight{kind}) where {T,kind} = ChebyshevWeight{kind,T}()
+AbstractQuasiVector{T}(::ChebyshevWeight{kind}) where {T,kind} = ChebyshevWeight{kind,T}()
+
+
 getproperty(w::ChebyshevWeight{1,T}, ::Symbol) where T = -one(T)/2
 getproperty(w::ChebyshevWeight{2,T}, ::Symbol) where T = one(T)/2
 
@@ -29,6 +33,8 @@ const ChebyshevUWeight = ChebyshevWeight{2}
 const ChebyshevT = Chebyshev{1}
 const ChebyshevU = Chebyshev{2}
 
+show(io::IO, P::ChebyshevTWeight{Float64}) = summary(io, P)
+show(io::IO, P::ChebyshevUWeight{Float64}) = summary(io, P)
 summary(io::IO, ::ChebyshevTWeight{Float64}) = print(io, "ChebyshevTWeight()")
 summary(io::IO, ::ChebyshevUWeight{Float64}) = print(io, "ChebyshevUWeight()")
 
@@ -71,6 +77,9 @@ chebysevuweight(d::AbstractInterval{T}) where T = ChebyshevUWeight{float(T)}[aff
 ==(::Chebyshev, ::Legendre) = false
 ==(::Legendre, ::Chebyshev) = false
 
+show(io::IO, w::ChebyshevT{Float64}) = summary(io, w)
+show(io::IO, w::ChebyshevU{Float64}) = summary(io, w)
+
 summary(io::IO, w::ChebyshevT{Float64}) = print(io, "ChebyshevT()")
 summary(io::IO, w::ChebyshevU{Float64}) = print(io, "ChebyshevU()")
 
@@ -143,28 +152,30 @@ recurrencecoefficients(C::ChebyshevU) = (Fill(2,∞), Zeros{Int}(∞), Ones{Int}
 # Mass matrix
 ###
 
-massmatrix(::ChebyshevT{V}) where V = Diagonal([convert(V,π); Fill(convert(V,π)/2,∞)])
-massmatrix(::ChebyshevU{V}) where V = Diagonal(Fill(convert(V,π)/2,∞))
+weightedgrammatrix(::ChebyshevT{V}) where V = Diagonal([convert(V,π); Fill(convert(V,π)/2,∞)])
+weightedgrammatrix(::ChebyshevU{V}) where V = Diagonal(Fill(convert(V,π)/2,∞))
 
-@simplify *(A::QuasiAdjoint{<:Any,<:Weighted{<:Any,<:ChebyshevT}}, B::ChebyshevT) = massmatrix(ChebyshevT{promote_type(eltype(A),eltype(B))}())
-@simplify *(A::QuasiAdjoint{<:Any,<:Weighted{<:Any,<:ChebyshevU}}, B::ChebyshevU) = massmatrix(ChebyshevU{promote_type(eltype(A),eltype(B))}())
+@simplify *(A::QuasiAdjoint{<:Any,<:Weighted{<:Any,<:ChebyshevT}}, B::ChebyshevT) = weightedgrammatrix(ChebyshevT{promote_type(eltype(A),eltype(B))}())
+@simplify *(A::QuasiAdjoint{<:Any,<:Weighted{<:Any,<:ChebyshevU}}, B::ChebyshevU) = weightedgrammatrix(ChebyshevU{promote_type(eltype(A),eltype(B))}())
 
-@simplify function *(A::QuasiAdjoint{<:Any,<:ChebyshevT}, B::ChebyshevT)
-    T = promote_type(eltype(A), eltype(B))
+function grammatrix(A::ChebyshevT{T}) where T
     f = (k,j) -> isodd(j-k) ? zero(T) : -(((1 + (-1)^(j + k))*(-1 + j^2 + k^2))/(j^4 + (-1 + k^2)^2 - 2j^2*(1 + k^2)))
     BroadcastMatrix{T}(f, 0:∞, (0:∞)')
 end
 
+@simplify function *(A::QuasiAdjoint{<:Any,<:ChebyshevT}, B::ChebyshevT)
+    T = promote_type(eltype(A), eltype(B))
+    grammatrix(ChebyshevT{T}())
+end
+
 @simplify function *(A::QuasiAdjoint{<:Any,<:ChebyshevT}, B::ChebyshevU)
     T = promote_type(eltype(A), eltype(B))
     f = (k,j) -> isodd(j-k) ? zero(T) : ((one(T) + (-1)^(j + k))*(1 + j))/((1 + j - k)*(1 + j + k))
     BroadcastMatrix{T}(f, 0:∞, (0:∞)')
 end
 
 
-
-@simplify function *(A::QuasiAdjoint{<:Any,<:Weighted{<:Any,<:ChebyshevU}}, B::Weighted{<:Any,<:ChebyshevU})
-    T = promote_type(eltype(A), eltype(B))
+function grammatrix(A::Weighted{T,<:ChebyshevU}) where T
     f = (k,j) ->  isodd(j-k) ? zero(T) : -((2*(one(T) + (-1)^(j + k))*(1 + j)*(1 + k))/((-1 + j - k)*(1 + j - k)*(1 + j + k)*(3 + j + k)))
     BroadcastMatrix{T}(f, 0:∞, (0:∞)')
 end
@@ -183,15 +194,14 @@ end
 ##########
 
 # Ultraspherical(1)\(D*Chebyshev())
-@simplify function *(D::Derivative{<:Any,<:ChebyshevInterval}, S::ChebyshevT)
-    T = promote_type(eltype(D),eltype(S))
-    A = _BandedMatrix((zero(T):∞)', ℵ₀, -1,1)
-    ApplyQuasiMatrix(*, ChebyshevU{T}(), A)
+function diff(S::ChebyshevT{T}; dims=1) where T
+    D = _BandedMatrix((zero(T):∞)', ℵ₀, -1,1)
+    ApplyQuasiMatrix(*, ChebyshevU{T}(), D)
 end
 
-@simplify function *(D::Derivative{<:Any,<:ChebyshevInterval}, W::Weighted{<:Any,<:ChebyshevU})
-    T = promote_type(eltype(D),eltype(W))
-    Weighted(ChebyshevT{T}()) * _BandedMatrix((-one(T):-one(T):(-∞))', ℵ₀, 1,-1)
+function diff(W::Weighted{T,<:ChebyshevU}; dims=1) where T
+    D =  _BandedMatrix((-one(T):-one(T):(-∞))', ℵ₀, 1,-1)
+    ApplyQuasiMatrix(*, Weighted(ChebyshevT{T}()), D)
 end
 
 

src/classical/chebyshev.jl

@@ -147,14 +147,14 @@ end
     Diagonal(Fill(2convert(TV,π),(axes(B,2),)))
 end
 
-@simplify function *(D::Derivative, F::Fourier)
-    TV = promote_type(eltype(D),eltype(F))
-    Fourier{TV}()*_BlockArray(Diagonal(Vcat([reshape([0.0],1,1)], (1.0:∞) .* Fill([0 -one(TV); one(TV) 0], ∞))), (axes(F,2),axes(F,2)))
+function diff(F::Fourier{T}; dims=1) where T
+    D = _BlockArray(Diagonal(Vcat([reshape([zero(T)],1,1)], (one(T):∞) .* Fill([0 -one(T); one(T) 0], ∞))), (axes(F,2),axes(F,2)))
+    F * D
 end
 
-@simplify function *(D::Derivative, F::Laurent)
-    TV = promote_type(eltype(D),eltype(F))
-    Laurent{TV}() * Diagonal(PseudoBlockVector((((1:∞) .÷ 2) .* (1 .- 2 .* iseven.(1:∞))) * convert(TV,im), (axes(F,2),)))
+function diff(F::Laurent{T}; dims=1) where T
+    D = Diagonal(PseudoBlockVector((((one(real(T)):∞) .÷ 2) .* (1 .- 2 .* iseven.(1:∞))) * convert(T,im), (axes(F,2),)))
+    F * D
 end
 
 

src/classical/fourier.jl

@@ -6,6 +6,11 @@ is a quasi-vector representing `exp(-x^2)` on ℝ.
 struct HermiteWeight{T} <: Weight{T} end
 
 HermiteWeight() = HermiteWeight{Float64}()
+
+AbstractQuasiArray{T}(::HermiteWeight) where T = HermiteWeight{T}()
+AbstractQuasiVector{T}(::HermiteWeight) where T = HermiteWeight{T}()
+
+
 axes(::HermiteWeight{T}) where T = (Inclusion{T}(ℝ),)
 ==(::HermiteWeight, ::HermiteWeight) = true
 
@@ -25,6 +30,11 @@ broadcasted(::typeof(sqrt), H::HermiteWeight{T}) where T = H .^ (one(T)/2)
 
 struct Hermite{T} <: OrthogonalPolynomial{T} end
 Hermite() = Hermite{Float64}()
+
+AbstractQuasiArray{T}(::Hermite) where T = Hermite{T}()
+AbstractQuasiMatrix{T}(::Hermite) where T = Hermite{T}()
+
+
 orthogonalityweight(::Hermite{T}) where T = HermiteWeight{T}()
 
 ==(::Hermite, ::Hermite) = true
@@ -46,27 +56,17 @@ broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), ::HermiteWeight{T}, ::Hermite
 jacobimatrix(H::Hermite{T}) where T = Tridiagonal(Fill(one(T)/2,∞), Zeros{T}(∞), one(T):∞)
 recurrencecoefficients(H::Hermite{T}) where T = Fill{T}(2,∞), Zeros{T}(∞), zero(T):2:∞
 
-massmatrix(::Hermite{T}) where T = Diagonal(sqrt(convert(T,π)) .* convert(T,2) .^ (0:∞) .* gamma.(one(T):∞))
+weightedgrammatrix(::Hermite{T}) where T = Diagonal(sqrt(convert(T,π)) .* convert(T,2) .^ (0:∞) .* gamma.(one(T):∞))
 
-@simplify *(A::QuasiAdjoint{<:Any,<:Weighted{<:Any,<:Hermite}}, B::Hermite) = massmatrix(Hermite{promote_type(eltype(A),eltype(B))}())
+@simplify *(A::QuasiAdjoint{<:Any,<:Weighted{<:Any,<:Hermite}}, B::Hermite) = weightedgrammatrix(Hermite{promote_type(eltype(A),eltype(B))}())
 
 ##########
 # Derivatives
-##########
-
-@simplify function *(D::Derivative, H::Hermite)
-    T = promote_type(eltype(D),eltype(H))
-    D = _BandedMatrix((zero(T):2:∞)', ℵ₀, -1,1)
-    H*D
-end
-
-@simplify function *(D::Derivative, W::Weighted{<:Any,<:Hermite})
-    T = promote_type(eltype(D),eltype(W))
-    D = _BandedMatrix(Fill(-one(T), 1, ∞), ℵ₀, 1,-1)
-    W*D
-end
+#####
 
-@simplify function *(D::Derivative, Q::OrthonormalWeighted{<:Any,<:Hermite})
+diff(H::Hermite{T}; dims=1) where T = ApplyQuasiMatrix(*, H, _BandedMatrix((zero(T):2:∞)', ℵ₀, -1,1))
+diff(W::Weighted{T,<:Hermite}; dims=1) where T = ApplyQuasiMatrix(*, W, _BandedMatrix(Fill(-one(T), 1, ∞), ℵ₀, 1,-1))
+function diff(Q::OrthonormalWeighted{<:Any,<:Hermite}; dims=1)
     X = jacobimatrix(Q.P)
-    Q * Tridiagonal(-X.ev, X.dv, X.ev)
+    ApplyQuasiMatrix(*, Q, Tridiagonal(-X.ev, X.dv, X.ev))
 end