Weighted and HalfWeighted for simpler operators (#21)

* Add special support for Weighted OPs

* Use HalfWeighted to give some structure to derivative operators

* increase coverage

Author: dlfivefifty

src/ClassicalOrthogonalPolynomials.jl

@@ -149,6 +149,11 @@ recurrencecoefficients(P) = error("Override for $(typeof(P))")
 
 const WeightedOrthogonalPolynomial{T, A<:AbstractQuasiVector, B<:OrthogonalPolynomial} = WeightedBasis{T, A, B}
 
+function isorthogonalityweighted(wS::WeightedOrthogonalPolynomial)
+    w,S = wS.args
+    w == orthogonalityweight(S)
+end
+
 """
     singularities(f)
 

src/jacobi.jl

@@ -142,15 +142,46 @@ WeightedJacobi(a,b) = JacobiWeight(a,b) .* Jacobi(a,b)
 WeightedJacobi{T}(a,b) where T = JacobiWeight{T}(a,b) .* Jacobi{T}(a,b)
 
 
+"""
+    HalfWeighted{lr}(Jacobi(a,b))
+
+is equivalent to `JacobiWeight(a,0) .* Jacobi(a,b)` (`lr = :a`) or
+`JacobiWeight(0,b) .* Jacobi(a,b)` (`lr = :b`)
+"""
+struct HalfWeighted{lr, T, PP<:AbstractQuasiMatrix{T}} <: Basis{T}
+    P::PP
+end
+
+HalfWeighted{lr}(P) where lr = HalfWeighted{lr,eltype(P),typeof(P)}(P)
+
+axes(Q::HalfWeighted) = axes(Q.P)
+copy(Q::HalfWeighted) = Q
+
+==(A::HalfWeighted, B::HalfWeighted) = A.P == B.P
+
+convert(::Type{WeightedJacobi}, Q::HalfWeighted{:a,T}) where T = JacobiWeight(Q.P.a,zero(T)) .* Q.P
+convert(::Type{WeightedJacobi}, Q::HalfWeighted{:b,T}) where T = JacobiWeight(zero(T),Q.P.b) .* Q.P
+
+getindex(Q::HalfWeighted, x::Union{Number,AbstractVector}, jr::Union{Number,AbstractVector}) = convert(WeightedJacobi, Q)[x,jr]
+
+broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, Q::HalfWeighted) = Q * (Q.P \ (x .* Q.P))
+
+\(w_A::HalfWeighted, w_B::HalfWeighted) = convert(WeightedJacobi, w_A) \ convert(WeightedJacobi, w_B)
+\(w_A::HalfWeighted, B::AbstractQuasiArray) = convert(WeightedJacobi, w_A) \ B
+\(A::AbstractQuasiArray, w_B::HalfWeighted) = A \ convert(WeightedJacobi, w_B)
+
 axes(::AbstractJacobi{T}) where T = (Inclusion{T}(ChebyshevInterval{real(T)}()), oneto(∞))
 ==(P::Jacobi, Q::Jacobi) = P.a == Q.a && P.b == Q.b
 ==(P::Legendre, Q::Jacobi) = Jacobi(P) == Q
 ==(P::Jacobi, Q::Legendre) = P == Jacobi(Q)
 ==(A::WeightedJacobi, B::WeightedJacobi) = A.args == B.args
 ==(A::WeightedJacobi, B::Jacobi{T}) where T = A == JacobiWeight(zero(T),zero(T)).*B
 ==(A::WeightedJacobi, B::Legendre) = A == Jacobi(B)
-==(A::Jacobi{T}, B::WeightedJacobi) where T = JacobiWeight(zero(T),zero(T)).*A == B
 ==(A::Legendre, B::WeightedJacobi) = Jacobi(A) == B
+==(A::Jacobi{T}, B::WeightedJacobi) where T = JacobiWeight(zero(T),zero(T)).*A == B
+==(A::Legendre, B::Weighted{<:Any,<:AbstractJacobi}) = A == B.P
+==(A::Weighted{<:Any,<:AbstractJacobi}, B::Legendre) = A.P == B
+
 
 summary(io::IO, P::Jacobi) = print(io, "Jacobi($(P.a), $(P.b))")
 
@@ -290,6 +321,16 @@ function \(w_A::WeightedJacobi, B::Jacobi)
     w_A \ (JacobiWeight(zero(a),zero(b)) .* B)
 end
 
+function \(A::AbstractJacobi, w_B::WeightedJacobi)
+    Ã = Jacobi(A)
+    (A \ Ã) * (Ã \ w_B)
+end
+function \(w_A::WeightedJacobi, B::AbstractJacobi)
+    B̃ = Jacobi(B)
+    (w_A \ B̃) * (B̃ \ B)
+end
+
+
 function broadcastbasis(::typeof(+), w_A::WeightedJacobi, w_B::WeightedJacobi)
     wA,A = w_A.args
     wB,B = w_B.args
@@ -327,26 +368,46 @@ end
 ##########
 
 # Jacobi(a+1,b+1)\(D*Jacobi(a,b))
-@simplify function *(D::Derivative{<:Any,<:AbstractInterval}, S::Jacobi)
-    A = _BandedMatrix((((1:∞) .+ (S.a + S.b))/2)', ℵ₀, -1,1)
-    ApplyQuasiMatrix(*, Jacobi(S.a+1,S.b+1), A)
+@simplify *(D::Derivative{<:Any,<:AbstractInterval}, S::Jacobi) = Jacobi(S.a+1,S.b+1) * _BandedMatrix((((1:∞) .+ (S.a + S.b))/2)', ℵ₀, -1,1)
+
+@simplify function *(D::Derivative{<:Any,<:AbstractInterval}, WS::Weighted{<:Any,<:Jacobi})
+    # L_1^t
+    S = WS.P
+    a,b = S.a, S.b
+    if a == b == 0
+        D*S
+    else
+        Weighted(Jacobi(a-1, b-1)) * _BandedMatrix((-2*(1:∞))', ℵ₀, 1,-1)
+    end
 end
 
+#L_6^t
+@simplify function *(D::Derivative{<:Any,<:AbstractInterval}, WS::HalfWeighted{:a,<:Any,<:Jacobi})
+    S = WS.P
+    a,b = S.a, S.b
+    HalfWeighted{:a}(Jacobi(a-1,b+1)) * Diagonal(-(a:∞))
+end
+
+#L_6
+@simplify function *(D::Derivative{<:Any,<:AbstractInterval}, WS::HalfWeighted{:b,<:Any,<:Jacobi})
+    S = WS.P
+    a,b = S.a, S.b
+    HalfWeighted{:b}(Jacobi(a+1,b-1)) * Diagonal(b:∞)
+end
+
+
 # Jacobi(a-1,b-1)\ (D*w*Jacobi(a,b))
 @simplify function *(D::Derivative{<:Any,<:AbstractInterval}, WS::WeightedJacobi)
     w,S = WS.args
     a,b = S.a, S.b
-    if w.a == 0 && w.b == 0
+    if isorthogonalityweighted(WS) # L_1^t
+        D * Weighted(S)
+    elseif w.a == w.b == 0
         D*S
     elseif iszero(w.a) && w.b == b #L_6
-        A = Diagonal(b:∞)
-        ApplyQuasiMatrix(*, JacobiWeight(w.a,b-1) .* Jacobi(a+1,b-1), A)
+        D * HalfWeighted{:b}(S)
     elseif iszero(w.b) && w.a == a #L_6^t
-        A = Diagonal(-(a:∞))
-        ApplyQuasiMatrix(*, JacobiWeight(a-1,w.b) .* Jacobi(a-1,b+1), A)
-    elseif w.a == a && w.b == b # L_1^t
-        A = _BandedMatrix((-2*(1:∞))', ℵ₀, 1,-1)
-        ApplyQuasiMatrix(*, JacobiWeight(a-1,b-1) .* Jacobi(a-1, b-1), A)
+        D * HalfWeighted{:a}(S)
     elseif iszero(w.a)
         W = (JacobiWeight(w.a, b-1) .* Jacobi(a+1, b-1)) \ (D * (JacobiWeight(w.a,b) .* S))
         J = Jacobi(a+1,b) # range Jacobi

src/normalized.jl

@@ -146,6 +146,13 @@ show(io::IO, ::MIME"text/plain", Q::Normalized) = show(io, Q)
 
 
 
+
+"""
+    OrthonormalWeighted(P)
+
+is the orthonormal with respect to L^2 basis given by
+`sqrt.(orthogonalityweight(P)) .* Normalized(P)`.
+"""
 struct OrthonormalWeighted{T, PP<:AbstractQuasiMatrix{T}} <: Basis{T}
     P::Normalized{T, PP}
 end
@@ -164,4 +171,33 @@ function getindex(Q::OrthonormalWeighted, x::Union{Number,AbstractVector}, jr::U
     w = orthogonalityweight(Q.P)
     sqrt.(w[x]) .* Q.P[x,jr]
 end
-broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, Q::OrthonormalWeighted) = Q * (Q.P \ (x .* Q.P))
+broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, Q::OrthonormalWeighted) = Q * (Q.P \ (x .* Q.P))
+
+"""
+    Weighted(P)
+
+is equivalent to `orthogonalityweight(P) .* P`
+"""
+struct Weighted{T, PP<:AbstractQuasiMatrix{T}} <: Basis{T}
+    P::PP
+end
+
+axes(Q::Weighted) = axes(Q.P)
+copy(Q::Weighted) = Q
+
+==(A::Weighted, B::Weighted) = A.P == B.P
+
+convert(::Type{WeightedOrthogonalPolynomial}, P::Weighted) = orthogonalityweight(P.P) .* P.P
+
+function getindex(Q::Weighted, x::Union{Number,AbstractVector}, jr::Union{Number,AbstractVector})
+    w = orthogonalityweight(Q.P)
+    w[x] .* Q.P[x,jr]
+end
+broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, Q::Weighted) = Q * (Q.P \ (x .* Q.P))
+
+\(w_A::Weighted, w_B::Weighted) = convert(WeightedOrthogonalPolynomial, w_A) \ convert(WeightedOrthogonalPolynomial, w_B)
+\(w_A::Weighted, B::AbstractQuasiArray) = convert(WeightedOrthogonalPolynomial, w_A) \ B
+\(A::AbstractQuasiArray, w_B::Weighted) = A \ convert(WeightedOrthogonalPolynomial, w_B)
+
+@simplify *(Ac::QuasiAdjoint{<:Any,<:Weighted}, wB::Weighted) = 
+    convert(WeightedOrthogonalPolynomial, parent(Ac))' * convert(WeightedOrthogonalPolynomial, wB)

src/ultraspherical.jl

@@ -38,7 +38,7 @@ const WeightedUltraspherical{T} = WeightedBasis{T,<:UltrasphericalWeight,<:Ultra
 
 WeightedUltraspherical(λ) = UltrasphericalWeight(λ) .* Ultraspherical(λ)
 WeightedUltraspherical{T}(λ) where T = UltrasphericalWeight{T}(λ) .* Ultraspherical{T}(λ)
-
+orthogonalityweight(C::Ultraspherical) = UltrasphericalWeight(C.λ)
 
 ultrasphericalc(n::Integer, λ, z::Number) = Base.unsafe_getindex(Ultraspherical{promote_type(typeof(λ),typeof(z))}(λ), z, n+1)
 
@@ -47,6 +47,11 @@ ultrasphericalc(n::Integer, λ, z::Number) = Base.unsafe_getindex(Ultraspherical
 ==(::ChebyshevT, ::Ultraspherical) = false
 ==(C::Ultraspherical, ::ChebyshevU) = isone(C.λ)
 ==(::ChebyshevU, C::Ultraspherical) = isone(C.λ)
+==(P::Ultraspherical, Q::Jacobi) = isone(2P.λ) && Jacobi(P) == Q
+==(P::Jacobi, Q::Ultraspherical) = isone(2Q.λ) && P == Jacobi(Q)
+==(P::Ultraspherical, Q::Legendre) = isone(2P.λ)
+==(P::Legendre, Q::Ultraspherical) = isone(2Q.λ)
+
 
 ###
 # interrelationships
@@ -95,20 +100,30 @@ end
     ApplyQuasiMatrix(*, Ultraspherical{eltype(S)}(S.λ+1), A)
 end
 
-# Ultraspherical(λ-1)\ (D*w*Ultraspherical(λ))
-@simplify function *(D::Derivative{<:Any,<:AbstractInterval}, WS::WeightedUltraspherical)
-    w,S = WS.args
+# Ultraspherical(λ-1)\ (D*wUltraspherical(λ))
+@simplify function *(D::Derivative{<:Any,<:AbstractInterval}, WS::Weighted{<:Any,<:Ultraspherical})
+    S = WS.P
     λ = S.λ
     T = eltype(WS)
-    if iszero(w.λ)
-        D*S
-    elseif w.λ == λ == 1
+    if λ == 1
         A = _BandedMatrix((-(1:∞))', ℵ₀, 1,-1)
         ApplyQuasiMatrix(*, ChebyshevTWeight{T}() .* ChebyshevT{T}(), A)
-    elseif w.λ == λ
+    else
         n = (0:∞)
         A = _BandedMatrix((-one(T)/(2*(λ-1)) * ((n.+1) .* (n .+ (2λ-1))))', ℵ₀, 1,-1)
         ApplyQuasiMatrix(*, WeightedUltraspherical{T}(λ-1), A)
+    end
+end
+
+# Ultraspherical(λ-1)\ (D*w*Ultraspherical(λ))
+@simplify function *(D::Derivative{<:Any,<:AbstractInterval}, WS::WeightedUltraspherical)
+    w,S = WS.args
+    λ = S.λ
+    T = eltype(WS)
+    if iszero(w.λ)
+        D*S
+    elseif isorthogonalityweighted(WS) # weights match
+        D * Weighted(S)
     else
         error("Not implemented")
     end
@@ -180,7 +195,7 @@ function \(w_A::WeightedUltraspherical, w_B::WeightedUltraspherical)
 
     if wA == wB
         A \ B
-    elseif B.λ == A.λ+1 && wB.λ == wA.λ+1
+    elseif B.λ == A.λ+1 && wB.λ == wA.λ+1 # Lower
         λ = A.λ
         _BandedMatrix(Vcat(((2λ:∞) .* ((2λ+1):∞) ./ (4λ .* (λ+1:∞)))',
                             Zeros(1,∞),

test/test_jacobi.jl

@@ -1,5 +1,5 @@
-using ClassicalOrthogonalPolynomials, FillArrays, BandedMatrices, ContinuumArrays, QuasiArrays, LazyArrays, FastGaussQuadrature, Test
-import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMatrix
+using ClassicalOrthogonalPolynomials, FillArrays, BandedMatrices, ContinuumArrays, QuasiArrays, LazyArrays, LazyBandedMatrices, FastGaussQuadrature, Test
+import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMatrix, arguments, Weighted, HalfWeighted
 
 @testset "Jacobi" begin
     @testset "JacobiWeight" begin
@@ -251,7 +251,7 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
         U = ChebyshevU()
         JT = Jacobi(T)
         JU = Jacobi(U)
-        
+
         @testset "recurrence degenerecies" begin
             A,B,C = recurrencecoefficients(JT)
             @test A[1] == 0.5
@@ -343,7 +343,7 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
         @test norm((X*Min - Min*X')[1:n,1:n]) ≤ 1E-13
         β = X[n,n+1]*Mi[n+1,n+1]
         @test (x-y) * P[x,1:n]'Mi[1:n,1:n]*P[y,1:n] ≈ P[x,n:n+1]' * (X*Min - Min*X')[n:n+1,n:n+1] * P[y,n:n+1] ≈ P[x,n:n+1]' * [0 -β; β 0] * P[y,n:n+1]
-        
+
         @testset "extrapolation" begin
             x,y = 0.1,3.4
             @test (x-y) * P[x,1:n]'Mi[1:n,1:n]*Base.unsafe_getindex(P,y,1:n) ≈ P[x,n:n+1]' * [0 -β; β 0] * Base.unsafe_getindex(P,y,n:n+1)
@@ -353,4 +353,24 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
     @testset "special syntax" begin
         @test jacobip.(0:5, 0.1, 0.2, 0.3) == Jacobi(0.1, 0.2)[0.3, 1:6]
     end
+
+    @testset "Weighted/HalfWeighted" begin
+        x = axes(Legendre(),1)
+        D = Derivative(x)
+        a,b = 0.1,0.2
+        B = Jacobi(a,b)
+        A = Jacobi(a-1,b-1)
+        D_W = Weighted(A) \ (D * Weighted(B))
+        @test (A * (D_W * (B \ exp.(x))))[0.1] ≈ (-a*(1+0.1) + b*(1-0.1) + (1-0.1^2)) *exp(0.1)
+
+        D_a = HalfWeighted{:a}(Jacobi(a-1,b+1)) \ (D * HalfWeighted{:a}(B))
+        D_b = HalfWeighted{:b}(Jacobi(a+1,b-1)) \ (D * HalfWeighted{:b}(B))
+        @test (Jacobi(a-1,b+1) * (D_a * (B \ exp.(x))))[0.1] ≈ (-a + 1-0.1) *exp(0.1)
+        @test (Jacobi(a+1,b-1) * (D_b * (B \ exp.(x))))[0.1] ≈ (b + 1+0.1) *exp(0.1)
+
+        @test HalfWeighted{:a}(B) \ HalfWeighted{:a}(B) isa Eye
+        @test HalfWeighted{:a}(B) \ (JacobiWeight(a,0) .* B) isa Eye
+
+        @test HalfWeighted{:a}(B) \ (x .* HalfWeighted{:a}(B)) isa LazyBandedMatrices.Tridiagonal
+    end
 end

test/test_odes.jl

@@ -4,7 +4,7 @@ using ClassicalOrthogonalPolynomials, ContinuumArrays, QuasiArrays, BandedMatric
 import QuasiArrays: MulQuasiMatrix
 import ClassicalOrthogonalPolynomials: oneto
 import ContinuumArrays: MappedBasisLayout, MappedWeightedBasisLayout
-import LazyArrays: arguments, ApplyMatrix
+import LazyArrays: arguments, ApplyMatrix, colsupport, MemoryLayout
 import SemiseparableMatrices: VcatAlmostBandedLayout
 
 @testset "ODEs" begin
@@ -96,7 +96,7 @@ import SemiseparableMatrices: VcatAlmostBandedLayout
         D = Derivative(axes(S,1))
         X = Diagonal(Inclusion(axes(S,1)))
 
-        @test_broken (Legendre() \ S)*(S\(w.*S))
+        @test ((Legendre() \ S)*(S\(w.*S)))[1:10,1:10] ≈ (Legendre() \ (w .* S))[1:10,1:10]
         @test (Ultraspherical(3/2)\(D^2*(w.*S)))[1:10,1:10] ≈ diagm(0 => -(2:2:20))
     end
 

test/test_ultraspherical.jl

@@ -101,4 +101,20 @@ using ClassicalOrthogonalPolynomials, BandedMatrices, LazyArrays, Test
     @testset "special syntax" begin
         @test ultrasphericalc.(0:5, 2, 0.3) == Ultraspherical(2)[0.3, 1:6]
     end
+
+    @testset "Ultraspherical vs Chebyshev and Jacobi" begin
+        @test Ultraspherical(1) == ChebyshevU()
+        @test ChebyshevU() == Ultraspherical(1)
+        @test Ultraspherical(0) ≠ ChebyshevT()
+        @test ChebyshevT() ≠ Ultraspherical(0)
+        @test Ultraspherical(1) ≠ Jacobi(1/2,1/2)
+        @test Jacobi(1/2,1/2) ≠ Ultraspherical(1)
+        @test Ultraspherical(1/2) == Jacobi(0,0)
+        @test Ultraspherical(1/2) == Legendre()
+        @test Jacobi(0,0) == Ultraspherical(1/2)
+        @test Legendre() == Ultraspherical(1/2)
+
+        @test Ultraspherical(1/2) \ (JacobiWeight(0,0) .* Jacobi(0,0)) isa Diagonal
+        @test (JacobiWeight(0,0) .* Jacobi(0,0)) \ Ultraspherical(1/2) isa Diagonal
+    end
 end