PartialDerivative -> diff (#198)

* PartialDerivative -> diff

* AbsLaplacianPower -> abslaplacian

* rect tests pass

* add coordinate, disk tests pass

* Fix some tests, support laplacian on rectangles

* Tests pass

* Update Project.toml

* Update test_triangle.jl

Project.toml

@@ -1,6 +1,6 @@
 name = "MultivariateOrthogonalPolynomials"
 uuid = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
-version = "0.8"
+version = "0.9"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -25,21 +25,21 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-ArrayLayouts = "1.0.9"
+ArrayLayouts = "1.11"
 BandedMatrices = "1"
 BlockArrays = "1.0"
 BlockBandedMatrices = "0.13"
-ClassicalOrthogonalPolynomials = "0.14.1"
-ContinuumArrays = "0.18"
+ClassicalOrthogonalPolynomials = "0.15"
+ContinuumArrays = "0.19"
 DomainSets = "0.7"
 FastTransforms = "0.17"
 FillArrays = "1.0"
-HarmonicOrthogonalPolynomials = "0.6.3"
+HarmonicOrthogonalPolynomials = "0.7"
 InfiniteArrays = "0.15"
 InfiniteLinearAlgebra = "0.9"
 LazyArrays = "2.3.1"
-LazyBandedMatrices = "0.11.1"
-QuasiArrays = "0.11"
+LazyBandedMatrices = "0.11.3"
+QuasiArrays = "0.12"
 RecurrenceRelationships = "0.2"
 SpecialFunctions = "1, 2"
 StaticArrays = "1"

README.md

@@ -12,7 +12,7 @@ julia> using MultivariateOrthogonalPolynomials, StaticArrays, LinearAlgebra
 julia> P = JacobiTriangle()
 JacobiTriangle(0, 0, 0)
 
-julia> x,y = first.(axes(P,1)), last.(axes(P,1));
+julia> x,y = coordinates(P);
 
 julia> u = P * (P \ (exp.(x) .* cos.(y))) # Expand in Triangle OPs
 JacobiTriangle(0, 0, 0) * [1.3365085377830084, 0.5687967596428205, -0.22812040274224554, 0.07733064070637755, 0.016169744493985644, -0.08714886622738759, 0.00338435674992512, 0.01220019521126353, -0.016867598915573725, 0.003930461395801074  …  ]

examples/diskheat.jl

@@ -3,8 +3,7 @@ pyplot() # pyplot supports disks
 
 Z = Zernike(1)
 W = Weighted(Z)
-xy = axes(W,1)
-x,y = first.(xy),last.(xy)
+x,y = coordinates(W)
 Δ = Z \ Laplacian(xy) * W
 S = Z \ W
 

examples/diskhelmholtz.jl

@@ -14,9 +14,8 @@ pyplot()
 
 Z = Zernike(1)
 W = Weighted(Z) # w*Z
-xy = axes(Z, 1);
-x, y = first.(xy), last.(xy);
-Δ = Z \ (Laplacian(xy) * W)
+x, y = coordinates(W)
+Δ = Z \ laplacian(W)
 S = Z \ W # identity
 
 

examples/multivariatelanczos.jl

@@ -3,8 +3,7 @@ using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, Test
 P = Legendre()
 P² = RectPolynomial(P,P)
 p₀ = expand(P², 𝐱 -> 1)
-𝐱 = axes(P²,1)
-x,y = first.(𝐱),last.(𝐱)
+x,y = coordinates(P²)
 w = P²/P²\ (x-y).^2
 
 w .* p₀

src/MultivariateOrthogonalPolynomials.jl

@@ -6,12 +6,12 @@ using ClassicalOrthogonalPolynomials, FastTransforms, BlockBandedMatrices, Block
       LazyArrays, SpecialFunctions, LinearAlgebra, BandedMatrices, LazyBandedMatrices, ArrayLayouts,
       HarmonicOrthogonalPolynomials, RecurrenceRelationships
 
-import Base: axes, in, ==, +, -, /, *, ^, \, copy, copyto!, OneTo, getindex, size, oneto, all, resize!, BroadcastStyle, similar, fill!, setindex!, convert, show, summary
+import Base: axes, in, ==, +, -, /, *, ^, \, copy, copyto!, OneTo, getindex, size, oneto, all, resize!, BroadcastStyle, similar, fill!, setindex!, convert, show, summary, diff
 import Base.Broadcast: Broadcasted, broadcasted, DefaultArrayStyle
 import DomainSets: boundary, EuclideanDomain
 
 import QuasiArrays: LazyQuasiMatrix, LazyQuasiArrayStyle, domain
-import ContinuumArrays: @simplify, Weight, weight, grid, plotgrid, TransformFactorization, ExpansionLayout, plotvalues, unweighted, plan_transform, checkpoints, transform_ldiv, AbstractBasisLayout, basis_axes, Inclusion, grammatrix, weaklaplacian, layout_broadcasted
+import ContinuumArrays: @simplify, Weight, weight, grid, plotgrid, TransformFactorization, ExpansionLayout, plotvalues, unweighted, plan_transform, checkpoints, transform_ldiv, AbstractBasisLayout, basis_axes, Inclusion, grammatrix, weaklaplacian, layout_broadcasted, laplacian, abslaplacian, laplacian_axis, abslaplacian_axis, diff_layout, operatororder, broadcastbasis
 
 import ArrayLayouts: MemoryLayout, sublayout, sub_materialize
 import BlockArrays: block, blockindex, BlockSlice, viewblock, blockcolsupport, AbstractBlockStyle, BlockStyle
@@ -23,17 +23,31 @@ import InfiniteArrays: InfiniteCardinal, OneToInf
 
 import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis, pad, recurrencecoefficients, clenshaw, weightedgrammatrix, Clenshaw
 import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, Plan,
-                                          PartialDerivative, AngularMomentum, BlockOneTo, BlockRange1, interlace,
+                                          AngularMomentum, angularmomentum, BlockOneTo, BlockRange1, interlace,
                                           MultivariateOPLayout, AbstractMultivariateOPLayout, MAX_PLOT_BLOCKS
 
 export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial,
        UnitTriangle, UnitDisk,
        JacobiTriangle, TriangleWeight, WeightedTriangle,
        DunklXuDisk, DunklXuDiskWeight, WeightedDunklXuDisk,
        Zernike, ZernikeWeight, zerniker, zernikez,
-       PartialDerivative, Laplacian, AbsLaplacianPower, AngularMomentum,
+       AngularMomentum,
        RadialCoordinate, Weighted, Block, jacobimatrix, KronPolynomial, RectPolynomial,
-       grammatrix, oneto
+       grammatrix, oneto, coordinates, Laplacian, AbsLaplacian, laplacian, abslaplacian, angularmomentum, weaklaplacian
+
+
+laplacian_axis(::Inclusion{<:SVector{2}}, A; dims...) = diff(A, (2,0); dims...) + diff(A, (0, 2); dims...)
+abslaplacian_axis(::Inclusion{<:SVector{2}}, A; dims...) = -(diff(A, (2,0); dims...) + diff(A, (0, 2); dims...))
+coordinates(P) = (first.(axes(P,1)), last.(axes(P,1)))
+
+function diff_layout(::AbstractBasisLayout, a, (k,j)::NTuple{2,Int}; dims...)
+    (k < 0 || j < 0) && throw(ArgumentError("order must be non-negative"))
+    k == j == 0 && return a
+    ((k,j) == (1,0) || (k,j) == (0,1)) && return diff(a, Val((k,j)); dims...)
+    k ≥ j && return diff(diff(a, (1,0)), (k-1,j))
+    diff(diff(a, (0,1)), (k,j-1))
+end
+
 
 include("ModalInterlace.jl")
 include("clenshawkron.jl")

src/disk.jl

@@ -262,16 +262,14 @@ plan_transform(Z::Zernike{T}, (N,)::Tuple{Block{1}}, dims=1) where T = ZernikeTr
 # Laplacian
 ###
 
-@simplify function *(Δ::Laplacian, WZ::Weighted{<:Any,<:Zernike})
+function laplacian(WZ::Weighted{T,<:Zernike}; dims...) where T
     @assert WZ.P.a == 0 && WZ.P.b == 1
-    T = eltype(eltype(WZ))
     WZ.P * ModalInterlace{T}(broadcast(k ->  Diagonal(-cumsum(k:8:∞)), 4:4:∞), (ℵ₀,ℵ₀), (0,0))
 end
 
-@simplify function *(Δ::Laplacian, Z::Zernike)
+function laplacian(Z::Zernike{T}; dims...) where T
     a,b = Z.a,Z.b
     @assert a == 0
-    T = promote_type(eltype(eltype(Δ)),eltype(Z)) # TODO: remove extra eltype
     D = Derivative(Inclusion(ChebyshevInterval{T}())) 
     Δs = BroadcastVector{AbstractMatrix{T}}((C,B,A) -> 4(HalfWeighted{:b}(C)\(D*HalfWeighted{:b}(B)))*(B\(D*A)), Normalized.(Jacobi.(b+2,a:∞)), Normalized.(Jacobi.(b+1,(a+1):∞)), Normalized.(Jacobi.(b,a:∞)))
     Δ = ModalInterlace(Δs, (ℵ₀,ℵ₀), (-2,2))
@@ -282,9 +280,9 @@ end
 # Fractional Laplacian
 ###
 
-function *(L::AbsLaplacianPower, WZ::Weighted{<:Any,<:Zernike{<:Any}})
-    @assert axes(L,1) == axes(WZ,1) && WZ.P.a == 0 && WZ.P.b == L.α
-    WZ.P * Diagonal(WeightedZernikeFractionalLaplacianDiag{typeof(L.α)}(L.α))
+function abslaplacian(WZ::Weighted{<:Any,<:Zernike}, α; dims...)
+    @assert WZ.P.a == 0 && WZ.P.b == α
+    WZ.P * Diagonal(WeightedZernikeFractionalLaplacianDiag{typeof(α)}(α))
 end
 
 # gives the entries for the (negative!) fractional Laplacian (-Δ)^(α) times (1-r^2)^α * Zernike(α)

src/rect.jl

@@ -52,18 +52,11 @@ function jacobimatrix(::Val{2}, P::RectPolynomial)
     Y = jacobimatrix(B)
     KronTrav(Y, Eye{eltype(Y)}(∞))
 end
-@simplify function *(Dx::PartialDerivative{1}, P::RectPolynomial)
-    A,B = P.args
-    U,M = (Derivative(axes(A,1))*A).args
-    # We want I ⊗ D² as A ⊗ B means B * X * A'
-    RectPolynomial(U,B) * KronTrav(Eye{eltype(M)}(∞), M)
+function diff(P::KronPolynomial{N}, order::NTuple{N,Int}; dims...) where N
+    diffs = diff.(P.args, order)
+    RectPolynomial(basis.(diffs)...) * KronTrav(coefficients.(diffs)...)
 end
 
-@simplify function *(Dx::PartialDerivative{2}, P::RectPolynomial)
-    A,B = P.args
-    U,M = (Derivative(axes(B,1))*B).args
-    RectPolynomial(A,U) * KronTrav(M, Eye{eltype(M)}(∞))
-end
 
 function weaklaplacian(P::RectPolynomial)
     A,B = P.args
@@ -188,3 +181,6 @@ function layout_broadcasted(::Tuple{ExpansionLayout{KronOPLayout{2}},KronOPLayou
 
     P * ClenshawKron(C, (recurrencecoefficients(A), recurrencecoefficients(B)), (jacobimatrix(T), jacobimatrix(U)))
 end
+
+
+broadcastbasis(::typeof(+), A::KronPolynomial, B::KronPolynomial) =  KronPolynomial(broadcastbasis.(+, A.args, B.args)...)

src/rectdisk.jl

@@ -41,7 +41,8 @@ const WeightedDunklXuDisk{T} = WeightedBasis{T,<:DunklXuDiskWeight,<:DunklXuDisk
 
 WeightedDunklXuDisk(β) = DunklXuDiskWeight(β) .* DunklXuDisk(β)
 
-@simplify function *(Dx::PartialDerivative{1}, P::DunklXuDisk)
+function diff(P::DunklXuDisk, ::Val{(1,0)}; dims=1)
+    @assert dims == 1
     β = P.β
     n = mortar(Fill.(oneto(∞),oneto(∞)))
     k = mortar(Base.OneTo.(oneto(∞)))
@@ -53,14 +54,15 @@ WeightedDunklXuDisk(β) = DunklXuDiskWeight(β) .* DunklXuDisk(β)
     DunklXuDisk(β+1) * _BandedBlockBandedMatrix(dat', axes(k,1), (-1,1), (0,2))
 end
 
-@simplify function *(Dy::PartialDerivative{2}, P::DunklXuDisk)
+function diff(P::DunklXuDisk{T}, ::Val{(0,1)}; dims=1) where T
+    @assert dims == 1
     β = P.β
     k = mortar(Base.OneTo.(oneto(∞)))
-    T = promote_type(eltype(Dy), eltype(P)) # avoid bug in convert
     DunklXuDisk(β+1) * _BandedBlockBandedMatrix(((k .+ T(2β)) ./ 2)', axes(k,1), (-1,1), (-1,1))
 end
 
-@simplify function *(Dx::PartialDerivative{1}, w_P::WeightedDunklXuDisk)
+function diff(w_P::WeightedDunklXuDisk, ::Val{(1,0)}; dims=1)
+    @assert dims == 1
     wP, P = w_P.args
     @assert P.β == wP.β
     β = P.β
@@ -74,11 +76,11 @@ end
     WeightedDunklXuDisk(β-1) * _BandedBlockBandedMatrix(dat', axes(k,1), (1,-1), (2,0))
 end
 
-@simplify function *(Dy::PartialDerivative{2}, w_P::WeightedDunklXuDisk)
+function diff(w_P::WeightedDunklXuDisk{T}, ::Val{(0,1)}; dims=1) where T
+    @assert dims == 1
     wP, P = w_P.args
     @assert P.β == wP.β
     k = mortar(Base.OneTo.(oneto(∞)))
-    T = promote_type(eltype(Dy), eltype(P)) # avoid bug in convert
     WeightedDunklXuDisk(P.β-1) * _BandedBlockBandedMatrix((T(-2).*k)', axes(k,1), (1,-1), (1,-1))
 end
 
@@ -178,7 +180,7 @@ function jacobimatrix(::Val{2}, P::DunklXuDisk)
     _BandedBlockBandedMatrix(dat', axes(k,1), (1,1), (1,1))
 end
 
-@simplify function *(A::AngularMomentum, P::DunklXuDisk)
+function angularmomentum(P::DunklXuDisk)
     β = P.β
     n = mortar(Fill.(oneto(∞),oneto(∞)))
     k = mortar(Base.OneTo.(oneto(∞)))

src/triangle.jl

@@ -147,7 +147,8 @@ summary(io::IO, P::TriangleWeight) = print(io, "x^$(P.a)*y^$(P.b)*(1-x-y)^$(P.c)
 
 orthogonalityweight(P::JacobiTriangle) = TriangleWeight(P.a, P.b, P.c)
 
-@simplify function *(Dx::PartialDerivative{1}, P::JacobiTriangle)
+function diff(P::JacobiTriangle{T}, ::Val{(1,0)}; dims=1) where T
+    @assert dims == 1
     a,b,c = P.a,P.b,P.c
     n = mortar(Fill.(oneto(∞),oneto(∞)))
     k = mortar(Base.OneTo.(oneto(∞)))
@@ -158,16 +159,17 @@ orthogonalityweight(P::JacobiTriangle) = TriangleWeight(P.a, P.b, P.c)
     JacobiTriangle(a+1,b,c+1) * _BandedBlockBandedMatrix(dat', axes(k,1), (-1,1), (0,1))
 end
 
-@simplify function *(Dy::PartialDerivative{2}, P::JacobiTriangle)
+function diff(P::JacobiTriangle{T}, ::Val{(0,1)}; dims=1) where T
+    @assert dims == 1
     a,b,c = P.a,P.b,P.c
     k = mortar(Base.OneTo.(oneto(∞)))
-    T = promote_type(eltype(Dy), eltype(P)) # avoid bug in convert
     JacobiTriangle(a,b+1,c+1) * _BandedBlockBandedMatrix((k .+ convert(T, b+c))', axes(k,1), (-1,1), (-1,1))
 end
 
 # TODO: The derivatives below only hold for a, b, c > 0.
-@simplify function *(Dx::PartialDerivative{1}, P::WeightedTriangle)
-    a, b, c = P.P.a, P.P.b, P.P.c
+function diff(w_P::WeightedTriangle{T}, ::Val{(1,0)}; dims=1) where T
+    @assert dims == 1
+    a, b, c = w_P.P.a, w_P.P.b, w_P.P.c
     n = mortar(Fill.(oneto(∞),oneto(∞)))
     k = mortar(Base.OneTo.(oneto(∞)))    
     scale = -(2k .+ (b + c - 1))
@@ -177,19 +179,13 @@ end
     WeightedTriangle(a-1, b, c-1) * _BandedBlockBandedMatrix(dat', axes(k, 1), (1, -1), (1, 0))
 end
 
-@simplify function *(Dy::PartialDerivative{2}, P::WeightedTriangle)
-    a, b, c = P.P.a, P.P.b, P.P.c 
+function diff(w_P::WeightedTriangle{T}, ::Val{(0,1)}; dims=1) where T
+    @assert dims == 1
+    a, b, c = w_P.P.a, w_P.P.b, w_P.P.c 
     k = mortar(Base.OneTo.(oneto(∞)))
-    T = promote_type(eltype(Dy), eltype(P)) # avoid bug in convert
     WeightedTriangle(a, b-1, c-1) * _BandedBlockBandedMatrix(-one(T) * k', axes(k, 1), (1, -1), (1, -1))
 end
 
-# @simplify function *(Δ::Laplacian, P)
-#     _lap_mul(P, eltype(axes(P,1)))
-# end
-
-
-
 function grammatrix(A::JacobiTriangle)
     @assert A == JacobiTriangle()
     n = mortar(Fill.(oneto(∞),oneto(∞)))