Don't dispatch on leftendpoint/rightendpoint in Evaluation (#239)

* Don't dispatch on leftendpoint/rightendpoint in Evaluation

* Specialize ldiffbc/rdiffbc for polynomial spaces

* Add bvp test

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunOrthogonalPolynomials"
 uuid = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
-version = "0.6.25"
+version = "0.6.26"
 
 [deps]
 ApproxFunBase = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"
@@ -20,7 +20,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-ApproxFunBase = "0.8.1"
+ApproxFunBase = "0.8.16"
 ApproxFunBaseTest = "0.1"
 Aqua = "0.5"
 BandedMatrices = "0.16, 0.17"

src/ApproxFunOrthogonalPolynomials.jl

@@ -42,7 +42,9 @@ import ApproxFunBase: Fun, SubSpace, WeightSpace, NoSpace, HeavisideSpace,
                     block, blockstart, blockstop, blocklengths, isblockbanded,
                     pointscompatible, affine_setdiff, complexroots,
                     ℓ⁰, recα, recβ, recγ, ℵ₀, ∞, RectDomain,
-                    assert_integer, supportsinplacetransform, ContinuousSpace
+                    assert_integer, supportsinplacetransform, ContinuousSpace, SpecialEvalPtType,
+                    isleftendpoint, isrightendpoint, evaluation_point, boundaryfn, ldiffbc, rdiffbc,
+                    LeftEndPoint, RightEndPoint
 
 import DomainSets: Domain, indomain, UnionDomain, FullSpace, Point,
             Interval, ChebyshevInterval, boundary, rightendpoint, leftendpoint,

src/Spaces/Chebyshev/ChebyshevOperators.jl

@@ -13,17 +13,7 @@ recγ(::Type{T},::Chebyshev,k) where {T} = one(T)/2   # one(T) ensures we get co
 
 ## Evaluation
 
-Evaluation(S::Chebyshev,x::typeof(leftendpoint),o::Integer) =
-    ConcreteEvaluation(S,x,o)
-Evaluation(S::Chebyshev,x::typeof(rightendpoint),o::Integer) =
-    ConcreteEvaluation(S,x,o)
-Evaluation(S::Chebyshev,x::Number,o::Integer) = ConcreteEvaluation(S,x,o)
-
-Evaluation(S::NormalizedPolynomialSpace{<:Chebyshev},x::typeof(leftendpoint),o::Integer) =
-    ConcreteEvaluation(S,x,o)
-Evaluation(S::NormalizedPolynomialSpace{<:Chebyshev},x::typeof(rightendpoint),o::Integer) =
-    ConcreteEvaluation(S,x,o)
-Evaluation(S::NormalizedPolynomialSpace{<:Chebyshev},x::Number,o::Integer) = ConcreteEvaluation(S,x,o)
+Evaluation(S::MaybeNormalized{<:Chebyshev},x::Number,o::Integer) = ConcreteEvaluation(S,x,o)
 
 function evaluatechebyshev(n::Integer,x::T) where T<:Number
     if n == 1
@@ -49,7 +39,21 @@ function forwardrecurrence(::Type{T},S::Chebyshev,r::AbstractUnitRange{<:Integer
     first(r) == 0 ? v : v[r .+ 1]
 end
 
-function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)},j::Integer) where {DD<:IntervalOrSegment,RR}
+function getindex(op::ConcreteEvaluation{<:Chebyshev{<:IntervalOrSegment},<:SpecialEvalPtType}, j::Integer)
+    _getindex_eval_endpoints(op, evaluation_point(op), j)
+end
+function getindex(op::ConcreteEvaluation{<:Chebyshev{<:IntervalOrSegment},<:SpecialEvalPtType}, j::AbstractUnitRange)
+    _getindex_eval_endpoints(op, evaluation_point(op), j)
+end
+
+function _getindex_eval_endpoints(op, x, j)
+    if isleftendpoint(x)
+        _getindex_eval_leftendpoint(op, x, j)
+    else
+        _getindex_eval_rightendpoint(op, x, j)
+    end
+end
+function _getindex_eval_leftendpoint(op::ConcreteEvaluation{<:Chebyshev{<:IntervalOrSegment}}, x, j::Integer)
     T=eltype(op)
     if op.order == 0
         ifelse(isodd(j),  # right rule
@@ -60,8 +64,7 @@ function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)
         op[j:j][1]
     end
 end
-
-function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(rightendpoint)},j::Integer) where {DD<:IntervalOrSegment,RR}
+function _getindex_eval_rightendpoint(op::ConcreteEvaluation{<:Chebyshev{<:IntervalOrSegment}}, x, j::Integer)
     T=eltype(op)
     if op.order == 0
         one(T)
@@ -70,8 +73,7 @@ function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(rightendpoint
         op[j:j][1]
     end
 end
-
-function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)},k::AbstractUnitRange) where {DD<:IntervalOrSegment,RR}
+function _getindex_eval_leftendpoint(op::ConcreteEvaluation{<:Chebyshev{<:IntervalOrSegment}}, x, k::AbstractUnitRange)
     Base.require_one_based_indexing(k)
     T=eltype(op)
     x = op.x
@@ -94,9 +96,7 @@ function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(leftendpoint)
 
     scal!(cst,ret)
 end
-
-
-function getindex(op::ConcreteEvaluation{<:Chebyshev{DD,RR},typeof(rightendpoint)},k::AbstractUnitRange) where {DD<:IntervalOrSegment,RR}
+function _getindex_eval_rightendpoint(op::ConcreteEvaluation{<:Chebyshev{<:IntervalOrSegment}}, x, k::AbstractUnitRange)
     Base.require_one_based_indexing(k)
     T=eltype(op)
     x = op.x

src/Spaces/Jacobi/JacobiOperators.jl

@@ -30,17 +30,10 @@ getindex(T::ConcreteDerivative{<:Jacobi}, k::Integer, j::Integer) =
 
 # Evaluation
 
-Evaluation(S::Jacobi,x::typeof(leftendpoint),o::Integer) =
-    ConcreteEvaluation(S,x,o)
-Evaluation(S::Jacobi,x::typeof(rightendpoint),o::Integer) =
-    ConcreteEvaluation(S,x,o)
-Evaluation(S::Jacobi,x::Number,o::Integer) = ConcreteEvaluation(S,x,o)
-
-Evaluation(S::NormalizedPolynomialSpace{<:Jacobi},x::typeof(leftendpoint),o::Integer) =
-    ConcreteEvaluation(S,x,o)
-Evaluation(S::NormalizedPolynomialSpace{<:Jacobi},x::typeof(rightendpoint),o::Integer) =
-    ConcreteEvaluation(S,x,o)
-Evaluation(S::NormalizedPolynomialSpace{<:Jacobi},x::Number,o::Integer) = ConcreteEvaluation(S,x,o)
+Evaluation(S::MaybeNormalized{<:Jacobi},x::Number,o::Integer) = ConcreteEvaluation(S,x,o)
+
+ldiffbc(d::MaybeNormalized{<:Union{Chebyshev, Ultraspherical, Jacobi}},k) = Evaluation(d,LeftEndPoint,k)
+rdiffbc(d::MaybeNormalized{<:Union{Chebyshev, Ultraspherical, Jacobi}},k) = Evaluation(d,RightEndPoint,k)
 
 ## Integral
 

src/Spaces/PolynomialSpace.jl

@@ -374,26 +374,12 @@ function forwardrecurrence(::Type{T},S::Space,r::AbstractRange,x::Number) where
 end
 
 getindex(op::ConcreteEvaluation{<:PolynomialSpace}, k::Integer) = op[k:k][1]
-# avoid ambiguity
-for OP in (:leftendpoint, :rightendpoint)
-    @eval getindex(op::ConcreteEvaluation{<:PolynomialSpace,typeof($OP)},k::Integer) =
-        op[k:k][1]
-end
-
-function getindex(op::ConcreteEvaluation{J,typeof(leftendpoint)},kr::AbstractRange) where J<:PolynomialSpace
-    _getindex_evaluation(op, leftendpoint(domain(op)), kr)
-end
 
-function getindex(op::ConcreteEvaluation{J,typeof(rightendpoint)},kr::AbstractRange) where J<:PolynomialSpace
-    _getindex_evaluation(op, rightendpoint(domain(op)), kr)
-end
-
-function getindex(op::ConcreteEvaluation{J,TT}, kr::AbstractRange) where {J<:PolynomialSpace,TT<:Number}
-    _getindex_evaluation(op, op.x, kr)
-end
+_evalpt(op, x::Number) = x
+_evalpt(op, x::SpecialEvalPtType) = boundaryfn(x)(domain(op))
 
-function _getindex_evaluation(op::ConcreteEvaluation{<:PolynomialSpace}, x, kr::AbstractRange)
-    _getindex_evaluation(eltype(op), op.space, op.order, x, kr)
+function getindex(op::ConcreteEvaluation{<:PolynomialSpace},kr::AbstractRange)
+    _getindex_evaluation(eltype(op), op.space, op.order, _evalpt(op, evaluation_point(op)), kr)
 end
 
 function _getindex_evaluation(::Type{T}, sp, order, x, kr::AbstractRange) where {T}

src/Spaces/Ultraspherical/DirichletSpace.jl

@@ -143,72 +143,88 @@ conversion_rule(b::ChebyshevDirichlet,a::Chebyshev)=b
 ## Evaluation Functional
 
 
-bandwidths(B::ConcreteEvaluation{ChebyshevDirichlet{1,0,D,R},typeof(leftendpoint)}) where {D,R} = 0,0
-bandwidths(B::ConcreteEvaluation{ChebyshevDirichlet{1,0,D,R},typeof(rightendpoint)}) where {D,R} = 0,ℵ₀
-bandwidths(B::ConcreteEvaluation{ChebyshevDirichlet{0,1,D,R},typeof(leftendpoint)}) where {D,R} = 0,ℵ₀
-bandwidths(B::ConcreteEvaluation{ChebyshevDirichlet{0,1,D,R},typeof(rightendpoint)}) where {D,R} = 0,0
-bandwidths(B::ConcreteEvaluation{ChebyshevDirichlet{1,1,D,R},typeof(leftendpoint)}) where {D,R} = 0,1
-bandwidths(B::ConcreteEvaluation{ChebyshevDirichlet{1,1,D,R},typeof(rightendpoint)}) where {D,R} = 0,1
+function bandwidths(B::ConcreteEvaluation{<:ChebyshevDirichlet{1,0},<:SpecialEvalPtType})
+    pt = evaluation_point(B)
+    if isleftendpoint(pt)
+        0,0
+    else
+        0,ℵ₀
+    end
+end
+function bandwidths(B::ConcreteEvaluation{<:ChebyshevDirichlet{0,1},<:SpecialEvalPtType})
+    pt = evaluation_point(B)
+    if isleftendpoint(pt)
+        0,ℵ₀
+    else
+        0,0
+    end
+end
+bandwidths(B::ConcreteEvaluation{<:ChebyshevDirichlet{1,1},<:SpecialEvalPtType}) = 0,1
 
-function getindex(B::ConcreteEvaluation{ChebyshevDirichlet{1,0,D,R},typeof(leftendpoint)},kr::AbstractRange) where {D,R}
-    d = domain(B)
+function getindex(B::ConcreteEvaluation{<:ChebyshevDirichlet,<:SpecialEvalPtType}, kr::AbstractRange)
+    _getindex_eval_endpoints(B, evaluation_point(B), kr)
+end
+
+function _getindex_default(B::ConcreteEvaluation{<:ChebyshevDirichlet}, x, kr)
+    S = Space(domain(B))
+    O = Evaluation(S,x,B.order) * Conversion(domainspace(B),S)
+    O[kr]
+end
 
+function _getindex_eval_leftendpoint(B::ConcreteEvaluation{<:ChebyshevDirichlet{1,0}}, x, kr::AbstractRange)
     if B.order == 0
         eltype(B)[k==1 ? 1.0 : 0.0 for k=kr]
     else
-        (Evaluation(d,B.x,B.order)*Conversion(domainspace(B)))[kr]
+        _getindex_default(B, x, kr)
     end
 end
 
-function getindex(B::ConcreteEvaluation{ChebyshevDirichlet{1,0,D,R},typeof(rightendpoint)},kr::AbstractRange) where {D,R}
-    d = domain(B)
-
+function _getindex_eval_rightendpoint(B::ConcreteEvaluation{<:ChebyshevDirichlet{1,0}}, x, kr::AbstractRange)
     if B.order == 0
         eltype(B)[k==1 ? 1.0 : 2.0 for k=kr]
     else
-        (Evaluation(d,B.x,B.order)*Conversion(domainspace(B)))[kr]
+        _getindex_default(B, x, kr)
     end
 end
 
-function getindex(B::ConcreteEvaluation{ChebyshevDirichlet{0,1,D,R},typeof(leftendpoint)},kr::AbstractRange) where {D,R}
+function _getindex_eval_leftendpoint(B::ConcreteEvaluation{<:ChebyshevDirichlet{0,1}}, x, kr::AbstractRange)
     S = Space(domain(B))
 
     if B.order == 0
         eltype(B)[k==1 ? 1.0 : -(-1)^k*2.0 for k=kr]
     else
-        (Evaluation(S,B.x,B.order)*Conversion(domainspace(B),S))[kr]
+        _getindex_default(B, x, kr)
     end
 end
 
-function getindex(B::ConcreteEvaluation{ChebyshevDirichlet{0,1,D,R},typeof(rightendpoint)},kr::AbstractRange) where {D,R}
+function _getindex_eval_rightendpoint(B::ConcreteEvaluation{<:ChebyshevDirichlet{0,1}}, x, kr::AbstractRange)
     S = Space(domain(B))
 
-
     if B.order == 0
         eltype(B)[k==1 ? 1.0 : 0.0 for k=kr]
     else
-        (Evaluation(S,B.x,B.order)*Conversion(domainspace(B),S))[kr]
+        _getindex_default(B, x, kr)
     end
 end
 
 
-function getindex(B::ConcreteEvaluation{ChebyshevDirichlet{1,1,D,R},typeof(leftendpoint)},kr::AbstractRange) where {D,R}
+function _getindex_eval_leftendpoint(B::ConcreteEvaluation{<:ChebyshevDirichlet{1,1}}, x, kr::AbstractRange)
     S = Space(domain(B))
 
     if B.order == 0
         eltype(B)[k==1 ? 1.0 : (k==2 ? -1.0 : 0.0) for k=kr]
     else
-        getindex(Evaluation(S,B.x,B.order)*Conversion(domainspace(B),S),kr)
+        _getindex_default(B, x, kr)
     end
 end
 
-function getindex(B::ConcreteEvaluation{ChebyshevDirichlet{1,1,D,R},typeof(rightendpoint)},kr::AbstractRange) where {D,R}
+function _getindex_eval_rightendpoint(B::ConcreteEvaluation{<:ChebyshevDirichlet{1,1}}, x, kr::AbstractRange)
     S = Space(domain(B))
 
     if B.order == 0
         eltype(B)[k==1||k==2 ? 1.0 : 0.0 for k=kr]
     else
-        getindex(Evaluation(S,B.x,B.order)*Conversion(domainspace(B),S),kr)
+        _getindex_default(B, x, kr)
     end
 end
 

src/Spaces/Ultraspherical/UltrasphericalOperators.jl

@@ -389,14 +389,4 @@ end
 ReverseOrientation(S::Ultraspherical) = ReverseOrientationWrapper(NegateEven(S,reverseorientation(S)))
 Reverse(S::Ultraspherical) = ReverseWrapper(NegateEven(S,S))
 
-Evaluation(S::Ultraspherical,x::typeof(leftendpoint),o::Integer) =
-    ConcreteEvaluation(S,x,o)
-Evaluation(S::Ultraspherical,x::typeof(rightendpoint),o::Integer) =
-    ConcreteEvaluation(S,x,o)
-Evaluation(S::Ultraspherical,x::Number,o::Integer) = ConcreteEvaluation(S,x,o)
-
-Evaluation(S::NormalizedPolynomialSpace{<:Ultraspherical},x::typeof(leftendpoint),o::Integer) =
-    ConcreteEvaluation(S,x,o)
-Evaluation(S::NormalizedPolynomialSpace{<:Ultraspherical},x::typeof(rightendpoint),o::Integer) =
-    ConcreteEvaluation(S,x,o)
-Evaluation(S::NormalizedPolynomialSpace{<:Ultraspherical},x::Number,o::Integer) = ConcreteEvaluation(S,x,o)
+Evaluation(S::MaybeNormalized{<:Ultraspherical},x::Number,o::Integer) = ConcreteEvaluation(S,x,o)

test/ChebyshevTest.jl

@@ -412,9 +412,15 @@ using ApproxFunOrthogonalPolynomials: forwardrecurrence
                     D3p = @inferred ev(sp, ep(d), 3)
                     @test Number(D3p * f) ≈ f'''(ep(d))
                 end
+                @test ldirichlet(sp) * f ≈ f(leftendpoint(domain(f)))
+                @test rdirichlet(sp) * f ≈ f(rightendpoint(domain(f)))
+                @test lneumann(sp) * f ≈ f'(leftendpoint(domain(f)))
+                @test rneumann(sp) * f ≈ f'(rightendpoint(domain(f)))
+                @inferred bvp(sp)[1]
             end
         end
 
+
         @testset "ChebyshevDirichlet" begin
             function Evaluation2(sp::ChebyshevDirichlet,x,ord)
                 S=Space(domain(sp))

test/JacobiTest.jl

@@ -635,6 +635,11 @@ using Static
                     D3p = @inferred ev(sp, ep(d), 3)
                     @test Number(D3p * f) ≈ f'''(ep(d))
                 end
+                @test ldirichlet(sp) * f ≈ f(leftendpoint(domain(f)))
+                @test rdirichlet(sp) * f ≈ f(rightendpoint(domain(f)))
+                @test lneumann(sp) * f ≈ f'(leftendpoint(domain(f)))
+                @test rneumann(sp) * f ≈ f'(rightendpoint(domain(f)))
+                @inferred bvp(sp)[1]
             end
         end
     end

test/UltrasphericalTest.jl

@@ -202,6 +202,11 @@ using Static
                     D3p = @inferred ev(sp, ep(d), 3)
                     @test Number(D3p * f) ≈ f'''(ep(d))
                 end
+                @test ldirichlet(sp) * f ≈ f(leftendpoint(domain(f)))
+                @test rdirichlet(sp) * f ≈ f(rightendpoint(domain(f)))
+                @test lneumann(sp) * f ≈ f'(leftendpoint(domain(f)))
+                @test rneumann(sp) * f ≈ f'(rightendpoint(domain(f)))
+                @inferred bvp(sp)[1]
             end
         end
     end