PartialDerivative -> diff

Author: dlfivefifty

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"
@@ -29,17 +29,17 @@ ArrayLayouts = "1.0.9"
 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"
+QuasiArrays = "0.12"
 RecurrenceRelationships = "0.2"
 SpecialFunctions = "1, 2"
 StaticArrays = "1"

src/MultivariateOrthogonalPolynomials.jl

@@ -23,15 +23,15 @@ 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, 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,
+       Laplacian, AbsLaplacianPower, AngularMomentum,
        RadialCoordinate, Weighted, Block, jacobimatrix, KronPolynomial, RectPolynomial,
        grammatrix, oneto
 

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))

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

src/rectdisk.jl

@@ -41,8 +41,8 @@ const WeightedDunklXuDisk{T} = WeightedBasis{T,<:DunklXuDiskWeight,<:DunklXuDisk
 
 WeightedDunklXuDisk(β) = DunklXuDiskWeight(β) .* DunklXuDisk(β)
 
-@simplify function *(Dx::PartialDerivative{1}, P::DunklXuDisk)
-    β = P.β
+function diff(::DunklXuDisk, ::Val{(1,0)}; dims=1)
+    @assert dims == 1
     n = mortar(Fill.(oneto(∞),oneto(∞)))
     k = mortar(Base.OneTo.(oneto(∞)))
     dat = BlockBroadcastArray(hcat,
@@ -53,14 +53,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(::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 +75,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
 

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(w_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,15 +159,16 @@ 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(w_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)
+function diff(w_P::WeightedTriangle{T}, ::Val{(1,0)}; dims=1) where T
+    @assert dims == 1
     a, b, c = P.P.a, P.P.b, P.P.c
     n = mortar(Fill.(oneto(∞),oneto(∞)))
     k = mortar(Base.OneTo.(oneto(∞)))    
@@ -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)
+function diff(w_P::WeightedTriangle{T}, ::Val{(0,1)}; dims=1) where T
+    @assert dims == 1
     a, b, c = P.P.a, P.P.b, 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(∞)))

test/test_rect.jl

@@ -79,7 +79,7 @@ using Base: oneto
         U² = RectPolynomial(U, U)
         C² = RectPolynomial(C, C)
         𝐱 = axes(T²,1)
-        D_x,D_y = PartialDerivative{1}(𝐱),PartialDerivative{2}(𝐱)
+        D_x,D_y = Derivative(𝐱,(1,0)),Derivative(𝐱,(0,1))
         D_x*T²
         D_y*T²
         U²\D_x*T²
@@ -103,7 +103,7 @@ using Base: oneto
 
         @testset "strong form" begin
             𝐱 = axes(W²,1)
-            D_x,D_y = PartialDerivative{1}(𝐱),PartialDerivative{2}(𝐱)
+            D_x,D_y = Derivative(𝐱,(1,0)),Derivative{2}(𝐱,(0,1))
             Δ = Q²\(D_x^2 + D_y^2)*W²
 
             K = Block.(1:200); @time L = Δ[K,K]; @time qr(L);

test/test_rectdisk.jl

@@ -39,8 +39,8 @@ import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering, Ang
 
         @test (DunklXuDisk() \ WeightedDunklXuDisk(1.0))[Block.(1:N), Block.(1:N)] ≈ (WeightedDunklXuDisk(0.0) \ WeightedDunklXuDisk(1.0))[Block.(1:N), Block.(1:N)]
 
-        ∂x = PartialDerivative{1}(axes(P, 1))
-        ∂y = PartialDerivative{2}(axes(P, 1))
+        ∂x = Derivative(axes(P, 1), (1,0))
+        ∂y = Derivative(axes(P, 1), (0,1))
 
         Dx = Q \ (∂x * P)
         Dy = Q \ (∂y * P)
@@ -67,8 +67,8 @@ import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering, Ang
         @test A[Block.(1:N), Block.(1:N)] ≈ C
         @test A2[Block.(1:N), Block.(1:N)] ≈ (A^2)[Block.(1:N), Block.(1:N)] ≈ A[Block.(1:N), Block.(1:N)]^2
 
-        ∂x = PartialDerivative{1}(axes(WQ, 1))
-        ∂y = PartialDerivative{2}(axes(WQ, 1))
+        ∂x = Derivative(axes(WQ, 1), (1,0))
+        ∂y = Derivative(axes(WQ, 1), (0,1))
 
         wDx = WP \ (∂x * WQ)
         wDy = WP \ (∂y * WQ)

test/test_triangle.jl

@@ -135,8 +135,8 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR
         P = JacobiTriangle()
         𝐱 = axes(P,1)
 
-        ∂ˣ = PartialDerivative{1}(𝐱)
-        ∂ʸ = PartialDerivative{2}(𝐱)
+        ∂ˣ = Derivative(𝐱, (1,0))
+        ∂ʸ = Derivative(𝐱, (0,1))
 
         @test eltype(∂ˣ) == eltype(∂ʸ) == Float64
 
@@ -467,8 +467,8 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR
         P¹ = JacobiTriangle(1,1,1)
         𝐱 = axes(P,1)
         x,y = first.(𝐱),last.(𝐱)
-        ∂ˣ = PartialDerivative{1}(𝐱)
-        ∂ʸ = PartialDerivative{2}(𝐱)
+        ∂ˣ = Derivative(𝐱, (1,0))
+        ∂ʸ = Derivative(𝐱, (0,1))
         L1 = x .* ∂ʸ
         L2 = y .* ∂ˣ
         L = x .* ∂ʸ - y .* ∂ˣ
@@ -518,8 +518,8 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR
                     P = Weighted(JacobiTriangle(a, b, c))
                     Pf = expand(P, f)
                     𝐱 = axes(P, 1)
-                    ∂ˣ = PartialDerivative{1}(𝐱)
-                    ∂ʸ = PartialDerivative{2}(𝐱)
+                    ∂ˣ = Derivative(𝐱, (1,0))
+                    ∂ʸ = Derivative(𝐱, (0,1))
                     Pfx = ∂ˣ * Pf
                     Pfy = ∂ʸ * Pf
 
@@ -579,8 +579,8 @@ import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleR
         P = JacobiTriangle()
         W = Weighted(JacobiTriangle(1,1,1))
         𝐱 = axes(W,1)
-        ∂_x = PartialDerivative{1}(𝐱)
-        ∂_y = PartialDerivative{2}(𝐱)
+        ∂_x = Derivative(𝐱, (1,0))
+        ∂_y = Derivative(𝐱, (0,1))
         Δ = -((∂_x*W)'*(∂_x*W) + (∂_y*W)'*(∂_y*W))
         M = W'W
         f = expand(P, splat((x,y) -> exp(x*cos(y))))