MemoryLayout(Weighted) isa WeightedLayout

Author: dlfivefifty

src/ClassicalOrthogonalPolynomials.jl

@@ -45,7 +45,7 @@ export OrthogonalPolynomial, Normalized, orthonormalpolynomial, LanczosPolynomia
             WeightedUltraspherical, WeightedChebyshev, WeightedChebyshevT, WeightedChebyshevU, WeightedJacobi,
             ∞, Derivative, .., Inclusion, 
             chebyshevt, chebyshevu, legendre, jacobi,
-            legendrep, jacobip, ultrasphericalc, laguerrel,hermiteh,
+            legendrep, jacobip, ultrasphericalc, laguerrel,hermiteh, normalizedjacobip,
             jacobimatrix, jacobiweight, legendreweight, chebyshevtweight, chebyshevuweight
 
 if VERSION < v"1.6-"

src/jacobi.jl

@@ -132,6 +132,7 @@ computes the `n`-th Jacobi polynomial, orthogonal with
 respec to `(1-x)^a*(1+x)^b`, at `z`.
 """
 jacobip(n::Integer, a, b, z::Number) = Base.unsafe_getindex(Jacobi{promote_type(typeof(a), typeof(b), typeof(z))}(a,b), z, n+1)
+normalizedjacobip(n::Integer, a, b, z::Number) = Base.unsafe_getindex(Normalized(Jacobi{promote_type(typeof(a), typeof(b), typeof(z))}(a,b)), z, n+1)
 
 OrthogonalPolynomial(w::JacobiWeight) = Jacobi(w.a, w.b)
 orthogonalityweight(P::Jacobi) = JacobiWeight(P.a, P.b)
@@ -148,7 +149,7 @@ WeightedJacobi{T}(a,b) where T = JacobiWeight{T}(a,b) .* Jacobi{T}(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}
+struct HalfWeighted{lr, T, PP<:AbstractQuasiMatrix{T}} <: AbstractWeighted{T}
     P::PP
 end
 
@@ -159,16 +160,18 @@ 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
+convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:a,T,<:Jacobi}) where T = JacobiWeight(Q.P.a,zero(T)) .* Q.P
+convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:b,T,<:Jacobi}) where T = JacobiWeight(zero(T),Q.P.b) .* Q.P
+convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:a,T,<:Normalized}) where T = JacobiWeight(Q.P.P.a,zero(T)) .* Q.P
+convert(::Type{WeightedOrthogonalPolynomial}, Q::HalfWeighted{:b,T,<:Normalized}) where T = JacobiWeight(zero(T),Q.P.P.b) .* Q.P
 
-getindex(Q::HalfWeighted, x::Union{Number,AbstractVector}, jr::Union{Number,AbstractVector}) = convert(WeightedJacobi, Q)[x,jr]
+getindex(Q::HalfWeighted, x::Union{Number,AbstractVector}, jr::Union{Number,AbstractVector}) = convert(WeightedOrthogonalPolynomial, 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)
+\(w_A::HalfWeighted, w_B::HalfWeighted) = convert(WeightedOrthogonalPolynomial, w_A) \ convert(WeightedOrthogonalPolynomial, w_B)
+\(w_A::HalfWeighted, B::AbstractQuasiArray) = convert(WeightedOrthogonalPolynomial, w_A) \ B
+\(A::AbstractQuasiArray, w_B::HalfWeighted) = A \ convert(WeightedOrthogonalPolynomial, 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
@@ -370,16 +373,6 @@ end
 # Jacobi(a+1,b+1)\(D*Jacobi(a,b))
 @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})
@@ -395,6 +388,29 @@ end
     HalfWeighted{:b}(Jacobi(a+1,b-1)) * Diagonal(b:∞)
 end
 
+for ab in (:(:a), :(:b))
+    @eval @simplify function *(D::Derivative{<:Any,<:AbstractInterval}, WS::HalfWeighted{$ab,<:Any,<:Normalized})
+        P,M = arguments(ApplyLayout{typeof(*)}(), WS.P)
+        D * HalfWeighted{$ab}(P) * M
+    end
+end
+
+
+@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
+    elseif iszero(a)
+        D * HalfWeighted{:b}(S)
+    elseif iszero(b)
+        D * HalfWeighted{:a}(S)
+    else
+        Weighted(Jacobi(a-1, b-1)) * _BandedMatrix((-2*(1:∞))', ℵ₀, 1,-1)
+    end
+end
+
 
 # Jacobi(a-1,b-1)\ (D*w*Jacobi(a,b))
 @simplify function *(D::Derivative{<:Any,<:AbstractInterval}, WS::WeightedJacobi)

src/normalized.jl

@@ -173,12 +173,17 @@ function getindex(Q::OrthonormalWeighted, x::Union{Number,AbstractVector}, jr::U
 end
 broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, Q::OrthonormalWeighted) = Q * (Q.P \ (x .* Q.P))
 
+
+abstract type AbstractWeighted{T} <: Basis{T} end
+
+MemoryLayout(::Type{<:AbstractWeighted}) = WeightedBasisLayout()
+
 """
     Weighted(P)
 
 is equivalent to `orthogonalityweight(P) .* P`
 """
-struct Weighted{T, PP<:AbstractQuasiMatrix{T}} <: Basis{T}
+struct Weighted{T, PP<:AbstractQuasiMatrix{T}} <: AbstractWeighted{T}
     P::PP
 end
 
@@ -200,4 +205,6 @@ broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, Q::Weighted) =
 \(A::AbstractQuasiArray, w_B::Weighted) = A \ convert(WeightedOrthogonalPolynomial, w_B)
 
 @simplify *(Ac::QuasiAdjoint{<:Any,<:Weighted}, wB::Weighted) = 
-    convert(WeightedOrthogonalPolynomial, parent(Ac))' * convert(WeightedOrthogonalPolynomial, wB)
+    convert(WeightedOrthogonalPolynomial, parent(Ac))' * convert(WeightedOrthogonalPolynomial, wB)
+
+summary(io::IO, Q::Weighted) = print(io, "Weighted($(Q.P))")

test/test_chebyshev.jl

@@ -1,6 +1,6 @@
 using ClassicalOrthogonalPolynomials, QuasiArrays, ContinuumArrays, BandedMatrices, LazyArrays, 
         FastTransforms, ArrayLayouts, Test, FillArrays, Base64, BlockArrays, LazyBandedMatrices, ForwardDiff
-import ClassicalOrthogonalPolynomials: Clenshaw, recurrencecoefficients, clenshaw, paddeddata, jacobimatrix, oneto
+import ClassicalOrthogonalPolynomials: Clenshaw, recurrencecoefficients, clenshaw, paddeddata, jacobimatrix, oneto, Weighted
 import LazyArrays: ApplyStyle
 import QuasiArrays: MulQuasiMatrix
 import Base: OneTo
@@ -128,6 +128,12 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
         u = wT * (wT \ @.(exp(x)/sqrt(1-x^2)))
         @test u[0.1] ≈ exp(0.1)/sqrt(1-0.1^2)
 
+        @testset "Weighted" begin
+            WT = Weighted(ChebyshevT())
+            @test wT[0.1,1:10] ≈ WT[0.1,1:10]
+            @test WT \ (exp.(x) ./ sqrt.(1 .- x.^2)) ≈ wT \ (exp.(x) ./ sqrt.(1 .- x.^2))
+        end
+
         @testset "mapped" begin
             x = Inclusion(0..1)
             wT̃ = wT[2x .- 1, :]