Static numbers in jacobi (#172)

* static numbers in jacobi

* version bump to v0.6.3

* Add tests

* Fix import

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunOrthogonalPolynomials"
 uuid = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
-version = "0.6.2"
+version = "0.6.3"
 
 [deps]
 ApproxFunBase = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"
@@ -15,11 +15,12 @@ IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-ApproxFunBase = "0.7.56"
+ApproxFunBase = "0.7.58"
 ApproxFunBaseTest = "0.1"
 Aqua = "0.5"
 BandedMatrices = "0.16, 0.17"
@@ -33,6 +34,7 @@ IntervalSets = "0.5, 0.6, 0.7"
 LazyArrays = "0.22"
 Reexport = "0.2, 1"
 SpecialFunctions = "0.10, 1.0, 2"
+Static = "0.8"
 StaticArrays = "1"
 julia = "1.6"
 

src/ApproxFunOrthogonalPolynomials.jl

@@ -1,7 +1,8 @@
 module ApproxFunOrthogonalPolynomials
 using Base, LinearAlgebra, Reexport, BandedMatrices, BlockBandedMatrices,
             BlockArrays, FillArrays, FastTransforms, IntervalSets,
-            DomainSets, Statistics, SpecialFunctions, FastGaussQuadrature
+            DomainSets, Statistics, SpecialFunctions, FastGaussQuadrature,
+            Static
 
 @reexport using ApproxFunBase
 

src/Spaces/Jacobi/Jacobi.jl

@@ -13,9 +13,9 @@ struct Jacobi{D<:Domain,R,T} <: PolynomialSpace{D,R}
 end
 Jacobi(b::T,a::T,d::Domain) where {T} =
     Jacobi{typeof(d),promote_type(T,real(prectype(d)))}(b, a, d)
-Legendre(domain) = Jacobi(0,0,domain)
+Legendre(domain) = Jacobi(static(0),static(0),domain)
 Legendre() = Legendre(ChebyshevInterval())
-Jacobi(b,a,d::Domain) = Jacobi(promote(b,a)...,d)
+Jacobi(b,a,d::Domain) = Jacobi(promote(dynamic(b), dynamic(a))...,d)
 Jacobi(b,a,d) = Jacobi(b,a,Domain(d))
 Jacobi(b,a) = Jacobi(b,a,ChebyshevInterval())
 Jacobi(A::Ultraspherical) = Jacobi(order(A)-0.5,order(A)-0.5,domain(A))
@@ -24,7 +24,7 @@ Jacobi(A::Chebyshev) = Jacobi(-0.5,-0.5,domain(A))
 NormalizedJacobi(s...) = NormalizedPolynomialSpace(Jacobi(s...))
 NormalizedLegendre(d...) = NormalizedPolynomialSpace(Legendre(d...))
 
-normalization(::Type{T}, sp::Jacobi, k::Int) where T = FastTransforms.Anαβ(k, sp.a, sp.b)
+normalization(::Type{T}, sp::Jacobi, k::Int) where T = FastTransforms.Anαβ(k, dynamic(sp.a), dynamic(sp.b))
 
 function Ultraspherical(J::Jacobi)
     if J.a == J.b
@@ -48,12 +48,13 @@ spacescompatible(a::Jacobi,b::Jacobi) = a.a ≈ b.a && a.b ≈ b.b && domainscom
 
 isapproxinteger_addhalf(a) = isapproxinteger(a+0.5)
 isapproxinteger_addhalf(::Integer) = false
+isapproxinteger_addhalf(@nospecialize ::StaticInt) = false
 function canonicalspace(S::Jacobi)
     if isapproxinteger_addhalf(S.a) && isapproxinteger_addhalf(S.b)
         Chebyshev(domain(S))
     else
         # return space with parameters in (-1,0.]
-        Jacobi(mod(S.b,-1),mod(S.a,-1),domain(S))
+        Jacobi(mod(dynamic(S.b),-1),mod(dynamic(S.a),-1),domain(S))
     end
 end
 
@@ -120,8 +121,8 @@ jacobip(r::AbstractRange,α,β,x::Number) = jacobip(promote_type(typeof(α),type
 
 jacobip(::Type{T},n::Integer,α,β,v) where {T} = jacobip(T,n:n,α,β,v)[1]
 jacobip(n::Integer,α,β,v) = jacobip(n:n,α,β,v)[1]
-jacobip(::Type{T},n,S::Jacobi,v) where {T} = jacobip(T,n,S.a,S.b,v)
-jacobip(n,S::Jacobi,v) = jacobip(n,S.a,S.b,v)
+jacobip(::Type{T},n,S::Jacobi,v) where {T} = jacobip(T,n,dynamic(S.a),dynamic(S.b),v)
+jacobip(n,S::Jacobi,v) = jacobip(n,dynamic(S.a),dynamic(S.b),v)
 
 
 

src/Spaces/Jacobi/JacobiOperators.jl

@@ -18,10 +18,10 @@ else
 end
 
 
-rangespace(D::ConcreteDerivative{J}) where {J<:Jacobi}=Jacobi(D.space.b+D.order,D.space.a+D.order,domain(D))
-bandwidths(D::ConcreteDerivative{J}) where {J<:Jacobi}=-D.order,D.order
+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
 
-getindex(T::ConcreteDerivative{J},k::Integer,j::Integer) where {J<:Jacobi} =
+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))
 
 
@@ -57,10 +57,10 @@ function Integral(J::Jacobi,k::Number)
 end
 
 
-rangespace(D::ConcreteIntegral{J}) where {J<:Jacobi}=Jacobi(D.space.b-D.order,D.space.a-D.order,domain(D))
-bandwidths(D::ConcreteIntegral{J}) where {J<:Jacobi}=D.order,0
+rangespace(D::ConcreteIntegral{<:Jacobi}) = Jacobi(D.space.b-D.order,D.space.a-D.order,domain(D))
+bandwidths(D::ConcreteIntegral{<:Jacobi}) = D.order,0
 
-function getindex(T::ConcreteIntegral{J},k::Integer,j::Integer) where J<:Jacobi
+function getindex(T::ConcreteIntegral{<:Jacobi}, k::Integer, j::Integer)
     @assert T.order==1
     if k≥2 && j==k-1
         complexlength(domain(T))./(k+T.space.a+T.space.b-2)
@@ -80,7 +80,7 @@ function Volterra(S::Jacobi, order::Integer)
 end
 
 rangespace(V::ConcreteVolterra{J}) where {J<:Jacobi}=Jacobi(-1.0,0.0,domain(V))
-bandwidths(V::ConcreteVolterra{J}) where {J<:Jacobi}=1,0
+bandwidths(::ConcreteVolterra{<:Jacobi}) = 1,0
 
 function getindex(V::ConcreteVolterra{J},k::Integer,j::Integer) where J<:Jacobi
     d=domain(V)
@@ -141,32 +141,36 @@ function Conversion(L::Jacobi,M::Jacobi)
 
     if isapproxinteger(L.a-M.a) && isapproxinteger(L.b-M.b)
         dm=domain(M)
-        D=typeof(dm)
         if isapprox(M.a,L.a) && isapprox(M.b,L.b)
             ConversionWrapper(Operator(I,L))
-        elseif (isapprox(M.b,L.b+1) && isapprox(M.a,L.a)) || (isapprox(M.b,L.b) && isapprox(M.a,L.a+1))
+        elseif (isapprox(M.b,L.b+static(1)) && isapprox(M.a,L.a)) ||
+            (isapprox(M.b,L.b) && isapprox(M.a,L.a+static(1)))
             ConcreteConversion(L,M)
-        elseif M.b > L.b+1
-            ConversionWrapper(TimesOperator(Conversion(Jacobi(M.b-1,M.a,dm),M),Conversion(L,Jacobi(M.b-1,M.a,dm))))
+        elseif M.b > L.b+static(1)
+            ConversionWrapper(
+                TimesOperator(
+                    Conversion(Jacobi(M.b-static(1),M.a,dm),M),
+                    Conversion(L,Jacobi(M.b-static(1),M.a,dm))))
         else  #if M.a >= L.a+1
-            ConversionWrapper(TimesOperator(Conversion(Jacobi(M.b,M.a-1,dm),M),Conversion(L,Jacobi(M.b,M.a-1,dm))))
+            ConversionWrapper(
+                TimesOperator(
+                    Conversion(Jacobi(M.b,M.a-static(1),dm),M),
+                    Conversion(L,Jacobi(M.b,M.a-static(1),dm))))
         end
-    elseif L.a ≈ L.b ≈ 0. && M.a ≈ M.b ≈ 0.5
+    elseif L.a ≈ L.b ≈ 0 && M.a ≈ M.b ≈ 0.5
         Conversion(L,Ultraspherical(L),Ultraspherical(M),M)
-    elseif L.a ≈ L.b ≈ 0. && M.a ≈ M.b ≈ -0.5
+    elseif L.a ≈ L.b ≈ 0 && M.a ≈ M.b ≈ -0.5
         Conversion(L,Ultraspherical(L),Chebyshev(M),M)
-    elseif L.a ≈ L.b ≈ -0.5 && M.a ≈ M.b ≈ 0.5
-        Conversion(L,Chebyshev(L),Ultraspherical(M),M)
     else # L.a - M.a ≈ L.b - M.b
         error("Implement for $L → $M")
     end
 end
 
-bandwidths(C::ConcreteConversion{J1,J2}) where {J1<:Jacobi,J2<:Jacobi}=(0,1)
+bandwidths(::ConcreteConversion{<:Jacobi,<:Jacobi}) = (0,1)
 
 
 
-function Base.getindex(C::ConcreteConversion{J1,J2,T},k::Integer,j::Integer) where {J1<:Jacobi,J2<:Jacobi,T}
+function Base.getindex(C::ConcreteConversion{<:Jacobi,<:Jacobi,T},k::Integer,j::Integer) where {T}
     L=C.domainspace
     if L.b+1==C.rangespace.b
         if j==k
@@ -220,7 +224,8 @@ function Conversion(A::Jacobi,B::PolynomialSpace)
 end
 
 isequalminhalf(x) = x == -0.5
-isequalminhalf(::Integer) = false
+isequalminhalf(@nospecialize ::Integer) = false
+isequalminhalf(@nospecialize ::StaticInt) = false
 
 function Conversion(A::Jacobi,B::Chebyshev)
     if isequalminhalf(A.a) && isequalminhalf(A.b)
@@ -300,23 +305,23 @@ end
 
 
 
-bandwidths(C::ConcreteConversion{US,J}) where {US<:Chebyshev,J<:Jacobi} = 0,0
-bandwidths(C::ConcreteConversion{J,US}) where {US<:Chebyshev,J<:Jacobi} = 0,0
+bandwidths(::ConcreteConversion{<:Chebyshev,<:Jacobi}) = 0,0
+bandwidths(::ConcreteConversion{<:Jacobi,<:Chebyshev}) = 0,0
 
 
-bandwidths(C::ConcreteConversion{US,J}) where {US<:Ultraspherical,J<:Jacobi} = 0,0
-bandwidths(C::ConcreteConversion{J,US}) where {US<:Ultraspherical,J<:Jacobi} = 0,0
+bandwidths(::ConcreteConversion{<:Ultraspherical,<:Jacobi}) = 0,0
+bandwidths(::ConcreteConversion{<:Jacobi,<:Ultraspherical}) = 0,0
 
 #TODO: Figure out how to unify these definitions
-function getindex(C::ConcreteConversion{CC,J,T},k::Integer,j::Integer) where {J<:Jacobi,CC<:Chebyshev,T}
+function getindex(::ConcreteConversion{<:Chebyshev,<:Jacobi,T}, k::Integer, j::Integer) where {T}
     if j==k
         one(T)/jacobip(T,k-1,-one(T)/2,-one(T)/2,one(T))
     else
         zero(T)
     end
 end
 
-function BandedMatrix(S::SubOperator{T,ConcreteConversion{CC,J,T},Tuple{UnitRange{Int},UnitRange{Int}}}) where {J<:Jacobi,CC<:Chebyshev,T}
+function BandedMatrix(S::SubOperator{T,ConcreteConversion{CC,J,T},NTuple{2,UnitRange{Int}}}) where {J<:Jacobi,CC<:Chebyshev,T}
     ret=BandedMatrix(Zeros, S)
     kr,jr = parentindices(S)
     k=(kr ∩ jr)
@@ -328,15 +333,15 @@ function BandedMatrix(S::SubOperator{T,ConcreteConversion{CC,J,T},Tuple{UnitRang
 end
 
 
-function getindex(C::ConcreteConversion{J,CC,T},k::Integer,j::Integer) where {J<:Jacobi,CC<:Chebyshev,T}
+function getindex(::ConcreteConversion{<:Jacobi,<:Chebyshev,T}, k::Integer, j::Integer) where {T}
     if j==k
         jacobip(T,k-1,-one(T)/2,-one(T)/2,one(T))
     else
         zero(T)
     end
 end
 
-function BandedMatrix(S::SubOperator{T,ConcreteConversion{J,CC,T},Tuple{UnitRange{Int},UnitRange{Int}}}) where {J<:Jacobi,CC<:Chebyshev,T}
+function BandedMatrix(S::SubOperator{T,ConcreteConversion{J,CC,T},NTuple{2,UnitRange{Int}}}) where {J<:Jacobi,CC<:Chebyshev,T}
     ret=BandedMatrix(Zeros, S)
     kr,jr = parentindices(S)
     k=(kr ∩ jr)
@@ -348,7 +353,7 @@ function BandedMatrix(S::SubOperator{T,ConcreteConversion{J,CC,T},Tuple{UnitRang
 end
 
 
-function getindex(C::ConcreteConversion{US,J,T},k::Integer,j::Integer) where {US<:Ultraspherical,J<:Jacobi,T}
+function getindex(C::ConcreteConversion{<:Ultraspherical,<:Jacobi,T}, k::Integer, j::Integer) where {T}
     if j==k
         S=rangespace(C)
         jp=jacobip(T,k-1,S.a,S.b,one(T))
@@ -359,7 +364,7 @@ function getindex(C::ConcreteConversion{US,J,T},k::Integer,j::Integer) where {US
     end
 end
 
-function BandedMatrix(S::SubOperator{T,ConcreteConversion{US,J,T},Tuple{UnitRange{Int},UnitRange{Int}}}) where {US<:Ultraspherical,J<:Jacobi,T}
+function BandedMatrix(S::SubOperator{T,ConcreteConversion{US,J,T},NTuple{2,UnitRange{Int}}}) where {US<:Ultraspherical,J<:Jacobi,T}
     ret=BandedMatrix(Zeros, S)
     kr,jr = parentindices(S)
     k=(kr ∩ jr)
@@ -375,7 +380,7 @@ end
 
 
 
-function getindex(C::ConcreteConversion{J,US,T},k::Integer,j::Integer) where {US<:Ultraspherical,J<:Jacobi,T}
+function getindex(C::ConcreteConversion{<:Jacobi,<:Ultraspherical,T}, k::Integer, j::Integer) where {T}
     if j==k
         S=domainspace(C)
         jp=jacobip(T,k-1,S.a,S.b,one(T))
@@ -386,7 +391,7 @@ function getindex(C::ConcreteConversion{J,US,T},k::Integer,j::Integer) where {US
     end
 end
 
-function BandedMatrix(S::SubOperator{T,ConcreteConversion{J,US,T},Tuple{UnitRange{Int},UnitRange{Int}}}) where {US<:Ultraspherical,J<:Jacobi,T}
+function BandedMatrix(S::SubOperator{T,ConcreteConversion{J,US,T},NTuple{2,UnitRange{Int}}}) where {US<:Ultraspherical,J<:Jacobi,T}
     ret=BandedMatrix(Zeros, S)
     kr,jr = parentindices(S)
     k=(kr ∩ jr)
@@ -403,6 +408,7 @@ end
 
 isapproxminhalf(a) = a ≈ -0.5
 isapproxminhalf(::Integer) = false
+isapproxminhalf(@nospecialize ::StaticInt) = false
 
 function union_rule(A::Jacobi,B::Jacobi)
     if domainscompatible(A,B)

src/Spaces/Ultraspherical/Ultraspherical.jl

@@ -1,4 +1,3 @@
-
 export Ultraspherical, NormalizedUltraspherical
 
 #Ultraspherical Spaces
@@ -106,17 +105,17 @@ ones(S::Ultraspherical) = ones(Float64, S)
 
 ## Fast evaluation
 
-function Base.first(f::Fun{Ultraspherical{Int,D,R}}) where {D,R}
+function Base.first(f::Fun{<:Ultraspherical{Int}})
     n = length(f.coefficients)
     n == 0 && return zero(cfstype(f))
     n == 1 && return first(f.coefficients)
     foldr(-,coefficients(f,Chebyshev))
 end
 
-Base.last(f::Fun{Ultraspherical{Int,D,R}}) where {D,R} = reduce(+,coefficients(f,Chebyshev))
+Base.last(f::Fun{<:Ultraspherical{Int}}) = reduce(+,coefficients(f,Chebyshev))
 
-Base.first(f::Fun{Ultraspherical{O,D,R}}) where {O,D,R} = f(leftendpoint(domain(f)))
-Base.last(f::Fun{Ultraspherical{O,D,R}}) where {O,D,R} = f(rightendpoint(domain(f)))
+Base.first(f::Fun{<:Ultraspherical}) = f(leftendpoint(domain(f)))
+Base.last(f::Fun{<:Ultraspherical}) = f(rightendpoint(domain(f)))
 
 function Fun(::typeof(identity), s::Ultraspherical)
     d = domain(s)

test/JacobiTest.jl

@@ -9,6 +9,7 @@ using ApproxFunBaseTest: testbandedbelowoperator, testbandedoperator, testspace,
                     testfunctional
 using ApproxFunOrthogonalPolynomials: jacobip
 using StaticArrays: SVector
+using Static
 
 @verbose @testset "Jacobi" begin
     @testset "Basic" begin
@@ -50,9 +51,9 @@ using StaticArrays: SVector
         @testset for d in [-1..1, 0..1]
             f = Fun(x->x^2, Chebyshev(d))
             C = space(f)
-            for J1 = Any[Jacobi(-0.5, -0.5, d), Legendre(d),
-                            Jacobi(0.5, 0.5, d), Jacobi(2.5, 1.5, d)]
-                for J in [J1, NormalizedPolynomialSpace(J1)]
+            for J1 in (Jacobi(-0.5, -0.5, d), Legendre(d),
+                            Jacobi(0.5, 0.5, d), Jacobi(2.5, 1.5, d))
+                for J in (J1, NormalizedPolynomialSpace(J1))
                     g = Fun(f, J)
                     if !any(isnan, coefficients(g))
                         @test Conversion(C, J) * f ≈ g
@@ -61,6 +62,64 @@ using StaticArrays: SVector
             end
         end
 
+        @testset "inference tests" begin
+            #= Note all cases are inferred as of now,
+            but as the situation eveolves in the future, more @inferred tests
+            may be added
+            There are also issues with static float promotion, because of which
+            we don't use the floating point orders directly in tests
+            See https://github.com/SciML/Static.jl/issues/97
+            =#
+            @testset "Jacobi" begin
+                CLL = @inferred Conversion(Legendre(), Legendre())
+                @test convert(Number, CLL) == 1
+                CNLNL = @inferred Conversion(NormalizedLegendre(), NormalizedLegendre())
+                @test convert(Number, CNLNL) == 1
+                CLNL = @inferred Conversion(Legendre(), NormalizedLegendre())
+                @test CLNL * Fun(Legendre()) ≈ Fun(NormalizedLegendre())
+                CNLL = @inferred Conversion(NormalizedLegendre(), Legendre())
+                @test CNLL * Fun(NormalizedLegendre()) ≈ Fun(Legendre())
+            end
+
+            @testset "Chebyshev" begin
+                CCL = @inferred Conversion(Chebyshev(), Legendre())
+                @test CCL * Fun(Chebyshev()) ≈ Fun(Legendre())
+                CLC = @inferred Conversion(Legendre(), Chebyshev())
+                @test CLC * Fun(Legendre()) ≈ Fun(Chebyshev())
+
+                @inferred Conversion(Chebyshev(), Jacobi(static(-0.5), static(-0.5)))
+                CCJmhalf = Conversion(Chebyshev(), Jacobi(-0.5, -0.5))
+                @test CCJmhalf * Fun(Chebyshev()) ≈ Fun(Jacobi(-0.5,-0.5))
+                @inferred Conversion(Jacobi(static(-0.5),static(-0.5)), Chebyshev())
+                CJmhalfC = Conversion(Jacobi(-0.5,-0.5), Chebyshev())
+                @test CJmhalfC * Fun(Jacobi(-0.5,-0.5)) ≈ Fun(Chebyshev())
+
+                @inferred Conversion(Chebyshev(), Jacobi(static(0.5), static(0.5)))
+                CCJmhalf = Conversion(Chebyshev(), Jacobi(0.5, 0.5))
+                @test CCJmhalf * Fun(Chebyshev()) ≈ Fun(Jacobi(0.5,0.5))
+
+                CCJ1 = Conversion(Chebyshev(), Jacobi(1,1))
+                @test CCJ1 * Fun(Chebyshev()) ≈ Fun(Jacobi(1,1))
+
+                CCJmix = Conversion(Chebyshev(), Jacobi(0.5,1.5))
+                @test CCJmix * Fun(Chebyshev()) ≈ Fun(Jacobi(0.5,1.5))
+            end
+
+            @testset "Ultraspherical" begin
+                CUL = @inferred Conversion(Ultraspherical(static(0.5)), Legendre())
+                @test CUL * Fun(Ultraspherical(0.5)) ≈ Fun(Legendre())
+                CLU = @inferred Conversion(Legendre(), Ultraspherical(static(0.5)))
+                @test CLU * Fun(Legendre()) ≈ Fun(Ultraspherical(0.5))
+
+                @inferred Conversion(Ultraspherical(static(0.5)), Jacobi(static(1),static(1)))
+                CU0J1 = Conversion(Ultraspherical(0.5), Jacobi(1,1))
+                @test CU0J1 * Fun(Ultraspherical(0.5)) ≈ Fun(Jacobi(1,1))
+                @inferred Conversion(Jacobi(static(1),static(1)), Ultraspherical(static(2.5)))
+                CJ1U2 = Conversion(Jacobi(1,1), Ultraspherical(2.5))
+                @test CJ1U2 * Fun(Jacobi(1,1)) ≈ Fun(Ultraspherical(2.5))
+            end
+        end
+
         @testset "conversion between spaces" begin
             for u in (Ultraspherical(1), Chebyshev())
                 @test NormalizedPolynomialSpace(Jacobi(u)) ==