Support weighted interface

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.3.0"
+version = "0.3.1"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -28,7 +28,7 @@ ArrayLayouts = "0.6.2"
 BandedMatrices = "0.16.5"
 BlockArrays = "0.14.1, 0.15"
 BlockBandedMatrices = "0.10"
-ContinuumArrays = "0.6"
+ContinuumArrays = "0.6.2"
 DomainSets = "0.4, 0.5"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3"

src/ClassicalOrthogonalPolynomials.jl

@@ -30,7 +30,7 @@ import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, Infini
 import ContinuumArrays: Basis, Weight, basis, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization,
     inbounds_getindex, grid, transform, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, Expansion,
     AffineQuasiVector, AffineMap, WeightLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
-    checkpoints
+    checkpoints, weight, unweightedbasis
 import FastTransforms: Λ, forwardrecurrence, forwardrecurrence!, _forwardrecurrence!, clenshaw, clenshaw!,
                         _forwardrecurrence_next, _clenshaw_next, check_clenshaw_recurrences, ChebyshevGrid, chebyshevpoints
 
@@ -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
 
@@ -157,18 +158,21 @@ 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
+==(A::HalfWeighted{lr}, B::HalfWeighted{lr}) where lr = A.P == B.P
+==(A::HalfWeighted, B::HalfWeighted) = false
 
-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 +374,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 +389,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,18 @@ 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()
+ContinuumArrays.unweightedbasis(wP::AbstractWeighted) = wP.P
+
 """
     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
 
@@ -187,17 +193,18 @@ copy(Q::Weighted) = Q
 
 ==(A::Weighted, B::Weighted) = A.P == B.P
 
-convert(::Type{WeightedOrthogonalPolynomial}, P::Weighted) = orthogonalityweight(P.P) .* P.P
+weight(wP::Weighted) = orthogonalityweight(wP.P)
 
-function getindex(Q::Weighted, x::Union{Number,AbstractVector}, jr::Union{Number,AbstractVector})
-    w = orthogonalityweight(Q.P)
-    w[x] .* Q.P[x,jr]
-end
+convert(::Type{WeightedOrthogonalPolynomial}, P::Weighted) = weight(P) .* unweightedbasis(P)
+
+getindex(Q::Weighted, x::Union{Number,AbstractVector}, jr::Union{Number,AbstractVector}) = weight(Q)[x] .* unweightedbasis(Q)[x,jr]
 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)
+    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,13 @@ 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))
+            @test WT[:,1:20] \ (exp.(x) ./ sqrt.(1 .- x.^2)) ≈ (WT \ (exp.(x) ./ sqrt.(1 .- x.^2)))[1:20]
+        end
+
         @testset "mapped" begin
             x = Inclusion(0..1)
             wT̃ = wT[2x .- 1, :]

test/test_jacobi.jl

@@ -1,5 +1,5 @@
 using ClassicalOrthogonalPolynomials, FillArrays, BandedMatrices, ContinuumArrays, QuasiArrays, LazyArrays, LazyBandedMatrices, FastGaussQuadrature, Test
-import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMatrix, arguments, Weighted, HalfWeighted
+import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMatrix, arguments, Weighted, HalfWeighted, WeightedOrthogonalPolynomial
 
 @testset "Jacobi" begin
     @testset "JacobiWeight" begin
@@ -352,6 +352,7 @@ 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]
+        @test normalizedjacobip.(0:5, 0.1, 0.2, 0.3) == Normalized(Jacobi(0.1, 0.2))[0.3, 1:6]
     end
 
     @testset "Weighted/HalfWeighted" begin
@@ -372,5 +373,19 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
         @test HalfWeighted{:a}(B) \ (JacobiWeight(a,0) .* B) isa Eye
 
         @test HalfWeighted{:a}(B) \ (x .* HalfWeighted{:a}(B)) isa LazyBandedMatrices.Tridiagonal
+
+        @test (D * HalfWeighted{:a}(Normalized(B)) * (Normalized(B) \ exp.(x)))[0.1] ≈ (-a + 1-0.1)*(1-0.1)^(a-1) *exp(0.1)
+        @test (D * HalfWeighted{:b}(Normalized(B)) * (Normalized(B) \ exp.(x)))[0.1] ≈ (b + 1+0.1) * (1+0.1)^(b-1)*exp(0.1)
+
+        @test (D * Weighted(Jacobi(0,0.1)))[0.1,1:10] ≈ (D * HalfWeighted{:b}(Jacobi(0,0.1)))[0.1,1:10]
+        @test (D * Weighted(Jacobi(0.1,0)))[0.1,1:10] ≈ (D * HalfWeighted{:a}(Jacobi(0.1,0)))[0.1,1:10]
+
+        @test HalfWeighted{:a}(Jacobi(0.2,0.1)) ≠ HalfWeighted{:b}(Jacobi(0.2,0.1))
+        @test HalfWeighted{:a}(Jacobi(0.2,0.1)) == HalfWeighted{:a}(Jacobi(0.2,0.1))
+
+        @test convert(WeightedOrthogonalPolynomial, HalfWeighted{:a}(Normalized(Jacobi(0.1,0.2))))[0.1,1:10] ≈
+            HalfWeighted{:a}(Normalized(Jacobi(0.1,0.2)))[0.1,1:10]
+        @test convert(WeightedOrthogonalPolynomial, HalfWeighted{:b}(Normalized(Jacobi(0.1,0.2))))[0.1,1:10] ≈
+            HalfWeighted{:b}(Normalized(Jacobi(0.1,0.2)))[0.1,1:10]            
     end
 end