Real -> Number (#18)

Author: dlfivefifty

Project.toml

@@ -25,6 +25,7 @@ QuasiArrays = "0.0.2"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 
 [targets]
-test = ["Test"]
+test = ["Test","ForwardDiff"]

src/ContinuumArrays.jl

@@ -31,7 +31,7 @@ cardinality(::AbstractInterval) = ℵ₁
 *(ℵ::AlephInfinity) = ℵ
 
 
-checkindex(::Type{Bool}, inds::AbstractInterval, i::Real) = (leftendpoint(inds) <= i) & (i <= rightendpoint(inds))
+checkindex(::Type{Bool}, inds::AbstractInterval, i::Number) = (leftendpoint(inds) <= i) & (i <= rightendpoint(inds))
 checkindex(::Type{Bool}, inds::AbstractInterval, i::Inclusion) = i.domain ⊆ inds
 function checkindex(::Type{Bool}, inds::AbstractInterval, I::AbstractArray)
     @_inline_meta

src/bases/jacobi.jl

@@ -9,7 +9,7 @@ end
 JacobiWeight(b::T, a::V) where {T,V} = JacobiWeight{promote_type(T,V)}(b,a)
 
 axes(::AbstractJacobiWeight) = (Inclusion(ChebyshevInterval()),)
-function getindex(w::JacobiWeight, x::Real)
+function getindex(w::JacobiWeight, x::Number)
     x ∈ axes(w,1) || throw(BoundsError())
     (1-x)^w.a * (1+x)^w.b
 end

src/bases/orthogonalpolynomials.jl

@@ -2,25 +2,7 @@ abstract type OrthogonalPolynomial{T} <: Basis{T} end
 
 @inline jacobioperator(P::OrthogonalPolynomial) =
     materialize(applied(\, P, applied(*, Diagonal(axes(P,1)), P)))
-
-# function forwardrecurrence(::Type{T},S::Space,r::AbstractRange,x::Number) where T
-#     if isempty(r)
-#         return T[]
-#     end
-#     n=maximum(r)+1
-#     v=Vector{T}(undef, n)  # x may be complex
-#     if n > 0
-#         v[1]=1
-#         if n > 1
-#             v[2] = muladd(recA(T,S,0),x,recB(T,S,0))
-#             @inbounds for k=2:n-1
-#                 v[k+1]=muladd(muladd(recA(T,S,k-1),x,recB(T,S,k-1)),v[k],-recC(T,S,k-1)*v[k-1])
-#             end
-#         end
-#     end
-
-#     return v[r.+1]
-# end    
+  
 
 function forwardrecurrence!(v::AbstractVector{T}, b::AbstractVector, a::AbstractVector, c::AbstractVector, x) where T
     isempty(v) && return v
@@ -77,7 +59,7 @@ _vec(a) = vec(a)
 _vec(a::Adjoint{<:Any,<:AbstractVector}) = a'
 bands(J) = _vec.(J.data.args)
 
-function getindex(P::OrthogonalPolynomial{T}, x::Real, n::OneTo) where T
+function getindex(P::OrthogonalPolynomial{T}, x::Number, n::OneTo) where T
     J = jacobioperator(P)
     b,a,c = bands(J)
     forwardrecurrence!(similar(n,T),b,a,c,x)
@@ -93,10 +75,10 @@ function getindex(P::OrthogonalPolynomial{T}, x::AbstractVector, n::OneTo) where
     V
 end
 
-getindex(P::OrthogonalPolynomial, x::Real, n::AbstractVector{<:Integer}) =
+getindex(P::OrthogonalPolynomial, x::Number, n::AbstractVector{<:Integer}) =
     P[x,OneTo(maximum(n))][n]
 
 getindex(P::OrthogonalPolynomial, x::AbstractVector, n::AbstractVector{<:Integer}) =
     P[x,OneTo(maximum(n))][:,n]    
 
-getindex(P::OrthogonalPolynomial, x::Real, n::Real) = P[x,OneTo(n)][end]
+getindex(P::OrthogonalPolynomial, x::Number, n::Number) = P[x,OneTo(n)][end]

src/bases/splines.jl

@@ -12,7 +12,7 @@ axes(B::Spline{o}) where o =
     (Inclusion(first(B.points)..last(B.points)), OneTo(length(B.points)+o-1))
 ==(A::Spline{o}, B::Spline{o}) where o = A.points == B.points
 
-function getindex(B::LinearSpline{T}, x::Real, k::Int) where T
+function getindex(B::LinearSpline{T}, x::Number, k::Int) where T
     x ∈ axes(B,1) && 1 ≤ k ≤ size(B,2)|| throw(BoundsError())
 
     p = B.points
@@ -25,7 +25,7 @@ function getindex(B::LinearSpline{T}, x::Real, k::Int) where T
     return (x-p[k+1])/(p[k]-p[k+1]) # x ≥ p[k]
 end
 
-function getindex(B::HeavisideSpline{T}, x::Real, k::Int) where T
+function getindex(B::HeavisideSpline{T}, x::Number, k::Int) where T
     x ∈ axes(B,1) && 1 ≤ k ≤ size(B,2)|| throw(BoundsError())
 
     p = B.points

src/bases/ultraspherical.jl

@@ -1,6 +1,6 @@
 struct ChebyshevWeight{T} <: AbstractJacobiWeight{T} end
 
-function getindex(w::ChebyshevWeight, x::Real)
+function getindex(w::ChebyshevWeight, x::Number)
     x ∈ axes(w,1) || throw(BoundsError())
     1/sqrt(1-x^2)
 end
@@ -9,7 +9,7 @@ struct UltrasphericalWeight{T,Λ} <: AbstractJacobiWeight{T}
     λ::Λ
 end
 
-function getindex(w::UltrasphericalWeight, x::Real)
+function getindex(w::UltrasphericalWeight, x::Number)
     x ∈ axes(w,1) || throw(BoundsError())
     (1-x^2)^(w.λ-one(w.λ)/2)
 end

src/operators.jl

@@ -127,7 +127,7 @@ IndexStyle(::Type{<:DiracDelta}) = IndexLinear()
 
 ==(a::DiracDelta, b::DiracDelta) = a.axis == b.axis && a.x == b.x
 
-function getindex(δ::DiracDelta{T}, x::Real) where T
+function getindex(δ::DiracDelta{T}, x::Number) where T
     x ∈ δ.axis || throw(BoundsError())
     x == δ.x ? inv(zero(T)) : zero(T)
 end

test/runtests.jl

@@ -1,7 +1,8 @@
-using ContinuumArrays, QuasiArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, Test, InfiniteArrays
+using ContinuumArrays, QuasiArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, Test, InfiniteArrays, ForwardDiff
 import ContinuumArrays: ℵ₁, materialize, Chebyshev, Ultraspherical, jacobioperator, SimplifyStyle
 import QuasiArrays: SubQuasiArray, MulQuasiMatrix, Vec, Inclusion, QuasiDiagonal, LazyQuasiArrayApplyStyle, LmaterializeApplyStyle
 import LazyArrays: MemoryLayout, ApplyStyle, Applied, colsupport
+import ForwardDiff: Dual
 
 
 @testset "Inclusion" begin
@@ -380,4 +381,18 @@ end
     cfs = C \ [1; zeros(n-2); 2] # Chebyshev coefficients
     u = P[:,1:n]*cfs  # interpret in basis
     @test u[0.1] ≈ (3cos(0.1)sec(1) + csc(1)sin(0.1))/2
+end
+
+@testset "Auto-diff" begin
+    U = Ultraspherical(1)
+    C = Ultraspherical(2)
+
+    f = x -> Chebyshev{eltype(x)}()[x,5]
+    @test ForwardDiff.derivative(f,0.1) ≈ 4*U[0.1,4]
+    f = x -> Chebyshev{eltype(x)}()[x,5][1]
+    @test ForwardDiff.gradient(f,[0.1]) ≈ [4*U[0.1,4]]
+    @test ForwardDiff.hessian(f,[0.1]) ≈ [8*C[0.1,3]]
+
+    f = x -> Chebyshev{eltype(x)}()[x,1:5]
+    @test ForwardDiff.derivative(f,0.1) ≈ [0;(1:4).*U[0.1,1:4]]
 end