Add type parameter to Jacobi (#158)

* Add type parameter to Jacobi

* Add tests

* fix multiplication bug

* Version bump to v0.6.0

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunOrthogonalPolynomials"
 uuid = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
-version = "0.5.22"
+version = "0.6.0"
 
 [deps]
 ApproxFunBase = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"

src/Spaces/Jacobi/Jacobi.jl

@@ -5,15 +5,15 @@ export Jacobi, NormalizedJacobi, Legendre, NormalizedLegendre
 `Jacobi(b,a)` represents the space spanned by Jacobi polynomials `P_k^{(a,b)}`,
 which are orthogonal with respect to the weight `(1+x)^β*(1-x)^α`
 """
-struct Jacobi{D<:Domain,R} <: PolynomialSpace{D,R}
-    b::R
-    a::R
+struct Jacobi{D<:Domain,R,T} <: PolynomialSpace{D,R}
+    b::T
+    a::T
     domain::D
-    Jacobi{D,R}(b,a,d) where {D,R} = new{D,R}(b,a,d)
+    Jacobi{D,R}(b::T,a::T,d::D) where {D,R,T} = new{D,R,T}(b,a,d)
 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(0,0,domain)
 Legendre() = Legendre(ChebyshevInterval())
 Jacobi(b,a,d::Domain) = Jacobi(promote(b,a)...,d)
 Jacobi(b,a,d) = Jacobi(b,a,Domain(d))
@@ -38,16 +38,18 @@ NormalizedUltraspherical(NS::NormalizedPolynomialSpace{<:Jacobi}) =
 NormalizedJacobi(NS::NormalizedPolynomialSpace{<:Union{Ultraspherical, Chebyshev}}) =
     NormalizedPolynomialSpace(Jacobi(NS.space))
 
-Base.promote_rule(::Type{Jacobi{D,R1}},::Type{Jacobi{D,R2}}) where {D,R1,R2} =
-    Jacobi{D,promote_type(R1,R2)}
-convert(::Type{Jacobi{D,R1}},J::Jacobi{D,R2}) where {D,R1,R2} =
-    Jacobi{D,R1}(J.b,J.a,J.domain)
+Base.promote_rule(::Type{Jacobi{D,R1,T1}},::Type{Jacobi{D,R2,T2}}) where {D,R1,R2,T1,T2} =
+    Jacobi{D,promote_type(R1,R2),promote_type(T1,T2)}
+convert(::Type{Jacobi{D,R1,T1}},J::Jacobi{D,R2,T2}) where {D,R1,R2,T1,T2} =
+    Jacobi{D,R1}(T1(J.b)::T1, T1(J.a)::T1, J.domain)::Jacobi{D,R1,T1}
 
 
 spacescompatible(a::Jacobi,b::Jacobi) = a.a ≈ b.a && a.b ≈ b.b && domainscompatible(a,b)
 
+isapproxinteger_addhalf(a) = isapproxinteger(a+0.5)
+isapproxinteger_addhalf(::Integer) = false
 function canonicalspace(S::Jacobi)
-    if isapproxinteger(S.a+0.5) && isapproxinteger(S.b+0.5)
+    if isapproxinteger_addhalf(S.a) && isapproxinteger_addhalf(S.b)
         Chebyshev(domain(S))
     else
         # return space with parameters in (-1,0.]

src/Spaces/Jacobi/JacobiOperators.jl

@@ -98,7 +98,7 @@ for (Func,Len,Sum) in ((:DefiniteIntegral,:complexlength,:sum),(:DefiniteLineInt
     @eval begin
         $Func(S::Jacobi{<:IntervalOrSegment}) = $ConcFunc(S)
 
-        function getindex(Σ::$ConcFunc{Jacobi{D,R},T},k::Integer) where {D<:IntervalOrSegment,R,T}
+        function getindex(Σ::$ConcFunc{<:Jacobi{<:IntervalOrSegment},T}, k::Integer) where {T}
             dsp = domainspace(Σ)
 
             if dsp.b == dsp.a == 0
@@ -109,7 +109,7 @@ for (Func,Len,Sum) in ((:DefiniteIntegral,:complexlength,:sum),(:DefiniteLineInt
             end
         end
 
-        function bandwidths(Σ::$ConcFunc{Jacobi{D,R}}) where {D<:IntervalOrSegment,R}
+        function bandwidths(Σ::$ConcFunc{<:Jacobi{<:IntervalOrSegment}})
             if domainspace(Σ).b == domainspace(Σ).a == 0
                 0,0  # first entry
             else
@@ -187,10 +187,10 @@ end
 
 # return the space that has banded Conversion to the other
 function conversion_rule(A::Jacobi,B::Jacobi)
-    if !isapproxinteger(A.a-B.a) || !isapproxinteger(A.b-B.b)
-        NoSpace()
-    else
+    if isapproxinteger(A.a-B.a) && isapproxinteger(A.b-B.b)
         Jacobi(min(A.b,B.b),min(A.a,B.a),domain(A))
+    else
+        NoSpace()
     end
 end
 
@@ -212,8 +212,11 @@ function Conversion(A::Jacobi,B::PolynomialSpace)
              ConversionWrapper(TimesOperator(Conversion(J,B),Conversion(A,J)))
 end
 
+isequalminhalf(x) = x == -0.5
+isequalminhalf(::Integer) = false
+
 function Conversion(A::Jacobi,B::Chebyshev)
-    if A.a == A.b == -0.5
+    if isequalminhalf(A.a) && isequalminhalf(A.b)
         ConcreteConversion(A,B)
     elseif A.a == A.b == 0
         ConversionWrapper(
@@ -230,7 +233,7 @@ function Conversion(A::Jacobi,B::Chebyshev)
 end
 
 function Conversion(A::Chebyshev,B::Jacobi)
-    if B.a == B.b == -0.5
+    if isequalminhalf(B.a) && isequalminhalf(B.b)
         ConcreteConversion(A,B)
     elseif B.a == B.b == 0
         ConversionWrapper(
@@ -248,10 +251,10 @@ end
 
 
 function Conversion(A::Jacobi,B::Ultraspherical)
-    if A.a == A.b == -0.5
+    if isequalminhalf(A.a) && isequalminhalf(A.b)
         ConversionWrapper(Conversion(Chebyshev(domain(A)),B)*
             ConcreteConversion(A,Chebyshev(domain(A))))
-    elseif A.a == A.b == order(B)-0.5
+    elseif isequalminhalf(A.a - order(B)) && isequalminhalf(A.b - order(B))
         ConcreteConversion(A,B)
     elseif A.a == A.b == 0
         ConversionWrapper(
@@ -268,10 +271,10 @@ function Conversion(A::Jacobi,B::Ultraspherical)
 end
 
 function Conversion(A::Ultraspherical,B::Jacobi)
-    if B.a == B.b == -0.5
+    if isequalminhalf(B.a) && isequalminhalf(B.b)
         ConversionWrapper(ConcreteConversion(Chebyshev(domain(A)),B)*
             Conversion(A,Chebyshev(domain(A))))
-    elseif B.a == B.b == order(A)-0.5
+    elseif isequalminhalf(B.a - order(A)) && isequalminhalf(B.b - order(A))
         ConcreteConversion(A,B)
     elseif B.a == B.b == 0
         ConversionWrapper(
@@ -391,9 +394,8 @@ function BandedMatrix(S::SubOperator{T,ConcreteConversion{J,US,T},Tuple{UnitRang
 end
 
 
-
-
-
+isapproxminhalf(a) = a ≈ -0.5
+isapproxminhalf(::Integer) = false
 
 function union_rule(A::Jacobi,B::Jacobi)
     if domainscompatible(A,B)
@@ -414,14 +416,14 @@ end
 
 function union_rule(A::Chebyshev,B::Jacobi)
     if domainscompatible(A, B)
-        if isapprox(B.a,-0.5) && isapprox(B.b,-0.5)
+        if isapproxminhalf(B.a) && isapproxminhalf(B.b)
             # the spaces are the same
             A
         else
             union(Jacobi(A),B)
         end
     else
-        if isapprox(B.a,-0.5) && isapprox(B.b,-0.5)
+        if isapproxminhalf(B.a) && isapproxminhalf(B.b)
             union(A, Chebyshev(domain(B)))
         else
             NoSpace()
@@ -431,14 +433,14 @@ end
 function union_rule(A::Ultraspherical,B::Jacobi)
     m=order(A)
     if domainscompatible(A, B)
-        if isapprox(B.a,m-0.5) && isapprox(B.b,m-0.5)
+        if isapproxminhalf(B.a-m) && isapproxminhalf(B.b-m)
             # the spaces are the same
             A
         else
             union(Jacobi(A),B)
         end
     else
-        if isapprox(B.a,m-0.5) && isapprox(B.b,m-0.5)
+        if isapproxminhalf(B.a-m) && isapproxminhalf(B.b-m)
             union(A, Ultraspherical(m, domain(B)))
         else
             NoSpace()
@@ -449,21 +451,21 @@ end
 for (OPrule,OP) in ((:conversion_rule,:conversion_type), (:maxspace_rule,:maxspace))
     @eval begin
         function $OPrule(A::Chebyshev,B::Jacobi)
-            if B.a ≈ -0.5 && B.b ≈ -0.5
+            if isapproxminhalf(B.a) && isapproxminhalf(B.b)
                 # the spaces are the same
                 A
-            elseif isapproxinteger(B.a+0.5) && isapproxinteger(B.b+0.5)
+            elseif isapproxinteger_addhalf(B.a) && isapproxinteger_addhalf(B.b)
                 $OP(Jacobi(A),B)
             else
                 NoSpace()
             end
         end
         function $OPrule(A::Ultraspherical,B::Jacobi)
             m = order(A)
-            if B.a ≈ m-0.5 && B.b ≈ m-0.5
+            if isapproxminhalf(B.a - m) && isapproxminhalf(B.b - m)
                 # the spaces are the same
                 A
-            elseif isapproxinteger(B.a+0.5) && isapproxinteger(B.b+0.5)
+            elseif isapproxinteger_addhalf(B.a) && isapproxinteger_addhalf(B.b)
                 $OP(Jacobi(A),B)
             else
                 NoSpace()

src/Spaces/Ultraspherical/UltrasphericalOperators.jl

@@ -134,8 +134,12 @@ function Conversion(A::Chebyshev,B::Ultraspherical)
     end
 end
 
+
+isequalhalf(x) = x == 0.5
+isequalhalf(::Integer) = false
+
 function Conversion(A::Ultraspherical,B::Chebyshev)
-    if order(A) == 1//2
+    if isequalhalf(order(A))
         ConcreteConversion(A,B)
     else
         error("Not implemented")

test/JacobiTest.jl

@@ -462,4 +462,15 @@ using StaticArrays: SVector
             end
         end
     end
+
+    @testset "type inference" begin
+        x = @inferred ApproxFunBase.maxspace(Jacobi(1,1), Jacobi(2,2))
+        @test x == Jacobi(2,2)
+
+        x = @inferred ApproxFunBase.union_rule(Chebyshev(), Jacobi(1,1))
+        @test x == Jacobi(-0.5, -0.5)
+
+        x = @inferred ApproxFunBase.conversion_rule(Jacobi(1,1), Jacobi(2,2))
+        @test x == Jacobi(1,1)
+    end
 end

test/UltrasphericalTest.jl

@@ -166,6 +166,13 @@ using ApproxFunOrthogonalPolynomials: jacobip
         end
     end
 
+    @testset "Multiplication" begin
+        M = Multiplication(Fun(), Ultraspherical(0.5))
+        f = Fun(Ultraspherical(0.5))
+        f2 = Fun(x->x^2, Ultraspherical(0.5))
+        @test M * f ≈ f2
+    end
+
     @testset "Integral" begin
         d = 0..1
         A = @inferred Integral(0..1)

test/showtest.jl

@@ -4,13 +4,15 @@ using ApproxFunBase
 	@test repr(NormalizedChebyshev()) == "NormalizedChebyshev()"
 	@test repr(Ultraspherical(1)) == "Ultraspherical(1)"
 	@test repr(Legendre()) == "Legendre()"
-	@test repr(Jacobi(1,2)) == "Jacobi(1.0,2.0)"
+	@test repr(Jacobi(1,2)) == "Jacobi(1,2)"
+	@test repr(Jacobi(1.0,2.0)) == "Jacobi(1.0,2.0)"
 
 	@test repr(Chebyshev(0..1)) == "Chebyshev(0..1)"
 	@test repr(NormalizedChebyshev(0..1)) == "NormalizedChebyshev(0..1)"
 	@test repr(Ultraspherical(1,0..1)) == "Ultraspherical(1,0..1)"
 	@test repr(Legendre(0..1)) == "Legendre(0..1)"
-	@test repr(Jacobi(1,2,0..1)) == "Jacobi(1.0,2.0,0..1)"
+	@test repr(Jacobi(1,2,0..1)) == "Jacobi(1,2,0..1)"
+	@test repr(Jacobi(1.0,2.0,0..1)) == "Jacobi(1.0,2.0,0..1)"
 
 	io = IOBuffer()
 	@testset "Derivative" begin