Concrete operators in normalized spaces (#330)

* Concrete operators in normalized spaces

* Concrete derivatives for Ultraspherical and Jacobi

* Add tests

Author: jishnub

src/Spaces/Chebyshev/Chebyshev.jl

@@ -26,7 +26,7 @@ const NormalizedChebyshev{D<:Domain,R} = NormalizedPolynomialSpace{Chebyshev{D,
 NormalizedChebyshev() = NormalizedPolynomialSpace(Chebyshev())
 NormalizedChebyshev(d) = NormalizedPolynomialSpace(Chebyshev(d))
 
-function Base.getproperty(S::Chebyshev{<:Any,<:Any},v::Symbol)
+function getproperty(S::Chebyshev, v::Symbol)
     if v==:b || v==:a
         -0.5
     else

src/Spaces/Chebyshev/ChebyshevOperators.jl

@@ -212,23 +212,24 @@ end
 
 ## Derivative
 
-function Derivative(sp::Chebyshev{<:IntervalOrSegment}, order::Number)
+function Derivative(sp::MaybeNormalized{<:Chebyshev{<:IntervalOrSegment}}, order::Number)
     assert_integer(order)
     ConcreteDerivative(sp,order)
 end
 
 
-rangespace(D::ConcreteDerivative{<:Chebyshev{<:IntervalOrSegment}}) =
+rangespace(D::ConcreteDerivative{<:MaybeNormalized{<:Chebyshev{<:IntervalOrSegment}}}) =
     Ultraspherical(D.order,domain(D))
 
-bandwidths(D::ConcreteDerivative{<:Chebyshev{<:IntervalOrSegment}}) = -D.order,D.order
-Base.stride(D::ConcreteDerivative{<:Chebyshev{<:IntervalOrSegment}}) = D.order
+bandwidths(D::ConcreteDerivative{<:MaybeNormalized{<:Chebyshev{<:IntervalOrSegment}}}) = -D.order,D.order
+Base.stride(D::ConcreteDerivative{<:MaybeNormalized{<:Chebyshev{<:IntervalOrSegment}}}) = D.order
 
-isdiag(D::ConcreteDerivative{<:Chebyshev{<:IntervalOrSegment}}) = false
+isdiag(D::ConcreteDerivative{<:MaybeNormalized{<:Chebyshev{<:IntervalOrSegment}}}) = false
 
-function getindex(D::ConcreteDerivative{<:Chebyshev{<:IntervalOrSegment},<:Any,T},k::Integer,j::Integer) where {T}
+function getindex(D::ConcreteDerivative{<:Chebyshev{<:IntervalOrSegment}}, k::Integer, j::Integer)
     m=D.order
     d=domain(D)
+    T = eltype(D)
 
     if j==k+m
         C=strictconvert(T,pochhammer(one(T),m-1)/2*(4/complexlength(d))^m)

src/Spaces/Jacobi/JacobiOperators.jl

@@ -1,7 +1,7 @@
 ## Derivative
 
 # specialize Derivative so that this is type-inferred even without constant propagation
-Derivative(J::Jacobi) = ConcreteDerivative(J,1)
+Derivative(J::MaybeNormalized{<:Jacobi}) = ConcreteDerivative(J,1)
 @inline function _Derivative(J::Jacobi, k::Number)
     assert_integer(k)
     if k==1
@@ -20,9 +20,9 @@ else
 end
 
 
-rangespace(D::ConcreteDerivative{<:Jacobi}) = Jacobi(D.space.b+D.order,D.space.a+D.order,domain(D))
-bandwidths(D::ConcreteDerivative{<:Jacobi}) = -D.order,D.order
-isdiag(D::ConcreteDerivative{<:Jacobi}) = false
+rangespace(D::ConcreteDerivative{<:MaybeNormalized{<:Jacobi}}) = Jacobi(D.space.b+D.order,D.space.a+D.order,domain(D))
+bandwidths(D::ConcreteDerivative{<:MaybeNormalized{<:Jacobi}}) = -D.order,D.order
+isdiag(D::ConcreteDerivative{<:MaybeNormalized{<:Jacobi}}) = false
 
 getindex(T::ConcreteDerivative{<:Jacobi}, k::Integer, j::Integer) =
     j==k+1 ? eltype(T)((k+1+T.space.a+T.space.b)/complexlength(domain(T))) : zero(eltype(T))

src/Spaces/PolynomialSpace.jl

@@ -430,6 +430,11 @@ end
 
 @containsconstants NormalizedPolynomialSpace
 
+@inline function getproperty(N::NormalizedPolynomialSpace, v::Symbol)
+    ((v == :a) || (v == :b)) && return getproperty(getfield(N, :space), v)
+    getfield(N, v)
+end
+
 domain(S::NormalizedPolynomialSpace) = domain(S.space)
 canonicalspace(S::NormalizedPolynomialSpace) = S.space
 setdomain(NS::NormalizedPolynomialSpace, d::Domain) = NormalizedPolynomialSpace(setdomain(canonicalspace(NS), d))
@@ -505,8 +510,8 @@ end
 
 # Conversion
 
-bandwidths(C::ConcreteConversion{NormalizedPolynomialSpace{S,D,R},S}) where {S,D,R} = (0, 0)
-bandwidths(C::ConcreteConversion{S,NormalizedPolynomialSpace{S,D,R}}) where {S,D,R} = (0, 0)
+bandwidths(::ConcreteConversion{<:NormalizedPolynomialSpace{S},S}) where {S} = (0, 0)
+bandwidths(::ConcreteConversion{S,<:NormalizedPolynomialSpace{S}}) where {S} = (0, 0)
 
 function getindex(C::ConcreteConversion{NormalizedPolynomialSpace{S,D,R},S,T},k::Integer,j::Integer) where {S,D<:IntervalOrSegment,R,T}
     if j==k
@@ -641,3 +646,12 @@ function evaluate(f::AbstractVector, S::NormalizedPolynomialSpace, x...)
     f_csp = mul_coefficients(C, f)
     evaluate(f_csp, csp, x...)
 end
+
+# Methods for concrete operators in normalized spaces
+function getindex(D::ConcreteDerivative{<:NormalizedPolynomialSpace}, k::Integer, j::Integer)
+    sp = domainspace(D)
+    csp = canonicalspace(sp)
+    Dcsp = ConcreteDerivative(csp, D.order)
+    C = Conversion_normalizedspace(csp, Val(:backward))
+    Dcsp[k, j] * C[j, j]
+end

src/Spaces/Ultraspherical/UltrasphericalOperators.jl

@@ -57,7 +57,7 @@ end
 #Derivative(d::IntervalOrSegment)=Derivative(1,d)
 
 
-function Derivative(sp::Ultraspherical{LT,DD}, m::Number) where {LT,DD<:IntervalOrSegment}
+function Derivative(sp::MaybeNormalized{<:Ultraspherical{<:Any,<:IntervalOrSegment}}, m::Number)
     assert_integer(m)
     ConcreteDerivative(sp,m)
 end
@@ -75,15 +75,15 @@ function Integral(sp::Ultraspherical{<:Any,<:IntervalOrSegment}, m::Number)
 end
 
 
-rangespace(D::ConcreteDerivative{<:Ultraspherical{LT,DD}}) where {LT,DD<:IntervalOrSegment} =
+rangespace(D::ConcreteDerivative{<:MaybeNormalized{<:Ultraspherical{<:Any,<:IntervalOrSegment}}}) =
     Ultraspherical(order(domainspace(D))+D.order,domain(D))
 
-bandwidths(D::ConcreteDerivative{<:Ultraspherical{LT,DD}}) where {LT,DD<:IntervalOrSegment} = -D.order,D.order
-bandwidths(D::ConcreteIntegral{<:Ultraspherical{LT,DD}}) where {LT,DD<:IntervalOrSegment} = D.order,-D.order
-Base.stride(D::ConcreteDerivative{<:Ultraspherical{LT,DD}}) where {LT,DD<:IntervalOrSegment} = D.order
+bandwidths(D::ConcreteDerivative{<:MaybeNormalized{<:Ultraspherical{<:Any,<:IntervalOrSegment}}}) = -D.order,D.order
+bandwidths(D::ConcreteIntegral{<:MaybeNormalized{<:Ultraspherical{<:Any,<:IntervalOrSegment}}}) = D.order,-D.order
+Base.stride(D::ConcreteDerivative{<:MaybeNormalized{<:Ultraspherical{<:Any,<:IntervalOrSegment}}}) = D.order
 
-isdiag(D::ConcreteDerivative{<:Ultraspherical{<:Any,<:IntervalOrSegment}}) = false
-isdiag(D::ConcreteIntegral{<:Ultraspherical{<:Any,<:IntervalOrSegment}}) = false
+isdiag(D::ConcreteDerivative{<:MaybeNormalized{<:Ultraspherical{<:Any,<:IntervalOrSegment}}}) = false
+isdiag(D::ConcreteIntegral{<:MaybeNormalized{<:Ultraspherical{<:Any,<:IntervalOrSegment}}}) = false
 
 function getindex(D::ConcreteDerivative{<:Ultraspherical{TT,DD},K,T},
                k::Integer,j::Integer) where {TT,DD<:IntervalOrSegment,K,T}

src/fastops.jl

@@ -72,8 +72,8 @@ end
 
 
 
-function BandedMatrix(S::SubOperator{T,ConcreteDerivative{Chebyshev{DD,RR},K,T},
-                                                     NTuple{2,UnitRange{Int}}}) where {T,K,DD,RR}
+function BandedMatrix(S::SubOperator{T,<:ConcreteDerivative{<:Chebyshev,<:Any,T},
+                                                     NTuple{2,UnitRange{Int}}}) where {T}
 
     n,m = size(S)
     ret = BandedMatrix{eltype(S)}(undef, (n,m), bandwidths(S))
@@ -98,8 +98,8 @@ function BandedMatrix(S::SubOperator{T,ConcreteDerivative{Chebyshev{DD,RR},K,T},
 end
 
 
-function BandedMatrix(S::SubOperator{T,ConcreteDerivative{Ultraspherical{LT,DD,RR},K,T},
-                                                  NTuple{2,UnitRange{Int}}}) where {T,K,DD,RR,LT}
+function BandedMatrix(S::SubOperator{T,<:ConcreteDerivative{<:Ultraspherical,T},
+                                                  NTuple{2,UnitRange{Int}}}) where {T}
     n,m = size(S)
     ret = BandedMatrix{eltype(S)}(undef, (n,m), bandwidths(S))
     kr,jr = parentindices(S)
@@ -116,3 +116,18 @@ function BandedMatrix(S::SubOperator{T,ConcreteDerivative{Ultraspherical{LT,DD,R
     ret
 end
 
+function BandedMatrix(S::SubOperator{T,<:ConcreteDerivative{<:NormalizedPolynomialSpace,<:Any,T},
+                                                     NTuple{2,UnitRange{Int}}}) where {T}
+
+    D = parent(S)
+    sp = domainspace(D)
+    csp = canonicalspace(sp)
+    ind1, ind2 = parentindices(S)
+    B = BandedMatrix(view(ConcreteDerivative(csp, D.order), ind1, ind2))
+    data = BandedMatrices.bandeddata(B)
+    C = ConcreteConversion(sp, csp)
+    for (ind, col) in enumerate(ind2)
+        @views data[:, ind] .*= C[col, col]
+    end
+    B
+end

test/ChebyshevTest.jl

@@ -311,6 +311,12 @@ include("testutils.jl")
             else
                 Derivative(s2)
             end
+            @test D1 isa ApproxFunBase.ConcreteDerivative
+            @test D2 isa ApproxFunBase.ConcreteDerivative
+            @test !isdiag(D1)
+            @test !isdiag(D2)
+            @test D1[1,2] ≈ D2[1,2]
+            @test D1[1:10, 1:10] ≈ D2[1:10, 1:10]
             f = x -> 3x^2 + 5x
             f1 = Fun(f, s1)
             f2 = Fun(f, s2)

test/JacobiTest.jl

@@ -404,6 +404,33 @@ include("testutils.jl")
         @test D * Fun(x->x^3, Legendre(0..1)) ≈ Fun(x->6x, Legendre(0..1))
     end
 
+    @testset "derivative in normalized space" begin
+        s1 = NormalizedLegendre(-1..1)
+        s2 = NormalizedLegendre()
+        @test s1 == s2
+        D1 = if VERSION >= v"1.8"
+            @inferred Derivative(s1)
+        else
+            Derivative(s1)
+        end
+        D2 = if VERSION >= v"1.8"
+            @inferred Derivative(s2)
+        else
+            Derivative(s2)
+        end
+        @test D1 isa ApproxFunBase.ConcreteDerivative
+        @test D2 isa ApproxFunBase.ConcreteDerivative
+        @test !isdiag(D1)
+        @test !isdiag(D2)
+        @test D1[1,2] ≈ D2[1,2]
+        @test D1[1:10, 1:10] ≈ D2[1:10, 1:10]
+        f = x -> 3x^2 + 5x
+        f1 = Fun(f, s1)
+        f2 = Fun(f, s2)
+        @test f1 == f2
+        @test D1 * f1 == D2 * f2
+    end
+
     @testset "identity Fun for interval domains" begin
         for d in [1..2, -1..1, 10..20], s in Any[Legendre(d), Jacobi(1, 2, d), Jacobi(1.2, 2.3, d)]
             f = Fun(identity, s)

test/UltrasphericalTest.jl

@@ -401,6 +401,34 @@ include("testutils.jl")
         @test M * f ≈ f2
     end
 
+    @testset "Derivative in a normalized space" begin
+        s1 = NormalizedUltraspherical(1,-1..1)
+        s2 = NormalizedUltraspherical(1)
+        @test s1 == s2
+        D1 = if VERSION >= v"1.8"
+            @inferred Derivative(s1)
+        else
+            Derivative(s1)
+        end
+        D2 = if VERSION >= v"1.8"
+            @inferred Derivative(s2)
+        else
+            Derivative(s2)
+        end
+        @test D1 isa ApproxFunBase.ConcreteDerivative
+        @test D2 isa ApproxFunBase.ConcreteDerivative
+        @test !isdiag(D1)
+        @test !isdiag(D2)
+        @test D1[1,2] ≈ D2[1,2]
+        @test D1[1:10, 1:10] ≈ D2[1:10, 1:10]
+        f = x -> 3x^2 + 5x
+        f1 = Fun(f, s1)
+        f2 = Fun(f, s2)
+        @test f1 == f2
+        df = x -> 6x + 5
+        @test D1 * f1 == D2 * f2 ≈ Fun(df, s2)
+    end
+
     @testset "Integral" begin
         d = 0..1
         A = @inferred Integral(0..1)