Simplify Jacobi conversion branches (#310)

* simplify Jacobi conversion branches

* Fix order of operators

* Tests, and improved inference in conversion to/from Chebyshev

* Add inference tests

* comment out unused branch

Author: jishnub

Project.toml

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

src/ApproxFunOrthogonalPolynomials.jl

@@ -57,7 +57,7 @@ import ApproxFunBase: Fun, SubSpace, WeightSpace, NoSpace, HeavisideSpace,
                     LeftEndPoint, RightEndPoint, normalizedspace, promotedomainspace,
                     bandmatrices_eigen, SymmetricEigensystem, SkewSymmetricEigensystem,
                     mean, # differs from Statistics.mean after https://github.com/JuliaApproximation/ApproxFunBase.jl/pull/506
-                    israggedbelow
+                    israggedbelow, bandwidthssum
 
 import DomainSets: Domain, indomain, UnionDomain, FullSpace, Point,
             Interval, ChebyshevInterval, boundary, rightendpoint, leftendpoint,

src/Spaces/Jacobi/JacobiOperators.jl

@@ -144,8 +144,15 @@ function _conversion_shiftordersbyone(L::Jacobi, M::Jacobi)
     # Conversion(J, M) = Conversion(Jacobi(M.b, L.a, dm), Jacobi(M.b, M.a, dm))
     CLJ = [ConcreteConversion(Jacobi(b-1,L.a,dm), Jacobi(b, L.a, dm)) for b in M.b:-1:L.b+1]
     CJM = [ConcreteConversion(Jacobi(M.b,a-1,dm), Jacobi(M.b, a, dm)) for a in M.a:-1:L.a+1]
-    C = [CJM; CLJ]
-    return ConversionWrapper(TimesOperator(C), L, M)
+    bw = bandwidthssum(bandwidths, CLJ) .+ bandwidthssum(bandwidths, CJM)
+    bbw = bandwidthssum(blockbandwidths, CLJ) .+ bandwidthssum(blockbandwidths, CJM)
+    sbbw = bandwidthssum(subblockbandwidths, CLJ) .+ bandwidthssum(subblockbandwidths, CJM)
+    op1 = isempty(CJM) ? CLJ[1] : CJM[1]
+    opend = isempty(CLJ) ? CJM[end] : CLJ[end]
+    ts = (size(op1, 1), size(opend, 2))
+    # we deliberately use Operator{T} to improve type-stability in Conversion
+    C = Operator{eltype(eltype(CLJ))}[CJM; CLJ]
+    return ConversionWrapper(TimesOperator(C, bw, ts, bbw, sbbw), L, M)
 end
 
 ## Conversion
@@ -162,28 +169,60 @@ function Conversion(L::Jacobi,M::Jacobi)
     dm=domain(M)
     dl=domain(L)
 
-    if isapproxinteger(L.a-M.a) && isapproxinteger(L.b-M.b) && M.b >= L.b && M.a >= L.a
-        if isapprox(L.a,M.a) && isapprox(L.b,M.b)
-            return ConversionWrapper(Operator(I,L))
-        elseif (isapprox(L.b+1,M.b) && isapprox(L.a,M.a)) ||
-                (isapprox(L.b,M.b) && isapprox(L.a+1,M.a))
+    La, Lb = L.a, L.b
+    Ma, Mb = M.a, M.b
+
+    if isapproxinteger(La-Ma) && isapproxinteger(Lb-Mb) && Mb >= Lb && Ma >= La
+        if isapprox(La,Ma) && isapprox(Lb,Mb)
+            return ConversionWrapper(Operator(I,L), L, M)
+        elseif (isapprox(Lb+1,Mb) && isapprox(La,Ma)) ||
+                (isapprox(Lb,Mb) && isapprox(La+1,Ma))
             return ConcreteConversion(L,M)
-        elseif L.a ≈ L.b && isapproxminhalf(L.a) && M.a ≈ M.b
-            return Conversion(L,Chebyshev(dl),Ultraspherical(M),M)
-        elseif L.a ≈ L.b && M.a ≈ M.b && isapproxminhalf(M.a)
-            return Conversion(L,Ultraspherical(L),Chebyshev(dm),M)
-        elseif L.a ≈ L.b && M.a ≈ M.b
-            return Conversion(L,Ultraspherical(L),Ultraspherical(M),M)
+        elseif La ≈ Lb && isapproxminhalf(La) && Ma ≈ Mb
+            C = Chebyshev(dl)
+            Conv_LC = ConcreteConversion(L, C)
+            U = Ultraspherical(M)
+            Conv_UM = ConcreteConversion(U, M)
+            return ConversionWrapper(Conv_UM * Conversion(C, U) * Conv_LC, L, M)
+            # return Conversion(L,Chebyshev(dl),Ultraspherical(M),M)
+        # elseif La ≈ Lb && Ma ≈ Mb && isapproxminhalf(Ma)
+        #     C = Chebyshev(dm)
+        #     Conv_CM = ConcreteConversion(C, M)
+        #     U = Ultraspherical(L)
+        #     Conv_LU = ConcreteConversion(L, U)
+        #     return ConversionWrapper(Conv_CM * Conversion(U, C) * Conv_LU, L, M)
+        #     # return Conversion(L,Ultraspherical(L),Chebyshev(dm),M)
+        elseif La ≈ Lb && Ma ≈ Mb && !isapproxminhalf(Ma)
+            UL = Ultraspherical(L)
+            Conv_LU = ConcreteConversion(L, UL)
+            UM = Ultraspherical(M)
+            Conv_UM = ConcreteConversion(UM, M)
+            return ConversionWrapper(Conv_UM * ultraconv_nonequal(UL, UM) * Conv_LU, L, M)
+            # return Conversion(L,Ultraspherical(L),Ultraspherical(M),M)
         else
             return _conversion_shiftordersbyone(L, M)
         end
-    elseif isapproxhalfoddinteger(L.a - M.a) && isapproxhalfoddinteger(L.b - M.b)
-        if L.a ≈ L.b && M.a ≈ M.b && isapproxminhalf(M.a)
-            return Conversion(L,Ultraspherical(L),Chebyshev(dm),M)
-        elseif L.a ≈ L.b && isapproxminhalf(L.a) && M.a ≈ M.b && M.a >= L.a
-            return Conversion(L,Chebyshev(dl),Ultraspherical(M),M)
-        elseif L.a ≈ L.b && M.a ≈ M.b && M.a >= L.a
-            return Conversion(L,Ultraspherical(L),Ultraspherical(M),M)
+    elseif isapproxhalfoddinteger(La - Ma) && isapproxhalfoddinteger(Lb - Mb)
+        if La ≈ Lb && Ma ≈ Mb && isapproxminhalf(Ma)
+            C = Chebyshev(dm)
+            Conv_CM = ConcreteConversion(C, M)
+            U = Ultraspherical(L)
+            Conv_LU = ConcreteConversion(L, U)
+            return ConversionWrapper(Conv_CM * Conversion(U, C) * Conv_LU, L, M)
+            # return Conversion(L,Ultraspherical(L),Chebyshev(dm),M)
+        elseif La ≈ Lb && isapproxminhalf(La) && Ma ≈ Mb && Ma >= La
+            C = Chebyshev(dl)
+            Conv_LC = ConcreteConversion(L, C)
+            U = Ultraspherical(M)
+            Conv_UM = ConcreteConversion(U, M)
+            return ConversionWrapper(Conv_UM * Conversion(C, U) * Conv_LC, L, M)
+            # return Conversion(L,Chebyshev(dl),Ultraspherical(M),M)
+        elseif La ≈ Lb && Ma ≈ Mb && Ma >= La
+            UL = Ultraspherical(L)
+            Conv_LU = ConcreteConversion(L, UL)
+            UM = Ultraspherical(M)
+            Conv_UM = ConcreteConversion(UM, M)
+            return ConversionWrapper(Conv_UM * Conversion(UL, UM) * Conv_LU, L, M)
         end
     end
 
@@ -264,7 +303,7 @@ function Conversion(A::Jacobi, B::Chebyshev)
         ConversionWrapper(SpaceOperator(ConcreteConversion(Ultraspherical(A), B), A, B))
     elseif A.a == A.b
         US = Ultraspherical(A)
-        ConversionWrapper(SpaceOperator(TimesOperator(Conversion(US,B), ConcreteConversion(A,US)), A, B))
+        ConversionWrapper(TimesOperator(Conversion(US,B), ConcreteConversion(A,US)), A, B)
     else
         J = Jacobi(B)
         ConcreteConversion(J,B)*Conversion(A,J)
@@ -279,7 +318,7 @@ function Conversion(A::Chebyshev, B::Jacobi)
         ConversionWrapper(SpaceOperator(ConcreteConversion(A, Ultraspherical(B)), A, B))
     elseif B.a == B.b
         US = Ultraspherical(B)
-        ConversionWrapper(SpaceOperator(TimesOperator(ConcreteConversion(US,B), Conversion(A,US)), A, B))
+        ConversionWrapper(TimesOperator(ConcreteConversion(US,B), Conversion(A,US)), A, B)
     else
         J = Jacobi(A)
         Conversion(J,B)*ConcreteConversion(A,J)
@@ -291,16 +330,16 @@ end
     @assert domain(A) == domain(B)
     if isequalminhalf(A.a) && isequalminhalf(A.b)
         C = Chebyshev(domain(A))
-        ConversionWrapper(SpaceOperator(
-            TimesOperator(Conversion(C,B), ConcreteConversion(A,C)), A, B))
+        ConversionWrapper(
+            TimesOperator(Conversion(C,B), ConcreteConversion(A,C)), A, B)
     elseif isequalminhalf(A.a - order(B)) && isequalminhalf(A.b - order(B))
         ConcreteConversion(A,B)
     elseif A.a == A.b == 0
         ConversionWrapper(SpaceOperator(Conversion(Ultraspherical(A), B), A, B))
     elseif A.a == A.b
         US = Ultraspherical(A)
-        ConversionWrapper(SpaceOperator(
-            TimesOperator(Conversion(US,B), ConcreteConversion(A,US)), A, B))
+        ConversionWrapper(
+            TimesOperator(Conversion(US,B), ConcreteConversion(A,US)), A, B)
     else
         J = Jacobi(B)
         ConcreteConversion(J,B)*Conversion(A,J)
@@ -311,16 +350,16 @@ end
     @assert domain(A) == domain(B)
     if isequalminhalf(B.a) && isequalminhalf(B.b)
         C = Chebyshev(domain(B))
-        ConversionWrapper(SpaceOperator(
-            TimesOperator(ConcreteConversion(C, B), Conversion(A, C)), A, B))
+        ConversionWrapper(
+            TimesOperator(ConcreteConversion(C, B), Conversion(A, C)), A, B)
     elseif isequalminhalf(B.a - order(A)) && isequalminhalf(B.b - order(A))
         ConcreteConversion(A,B)
     elseif B.a == B.b == 0
         ConversionWrapper(SpaceOperator(Conversion(A, Ultraspherical(B)), A, B))
     elseif B.a == B.b
         US = Ultraspherical(B)
-        ConversionWrapper(SpaceOperator(
-            TimesOperator(ConcreteConversion(US,B), Conversion(A,US)), A, B))
+        ConversionWrapper(
+            TimesOperator(ConcreteConversion(US,B), Conversion(A,US)), A, B)
     else
         J = Jacobi(A)
         Conversion(J,B)*ConcreteConversion(A,J)

src/Spaces/Ultraspherical/UltrasphericalOperators.jl

@@ -166,6 +166,26 @@ end
 
 maxspace_rule(A::Ultraspherical,B::Chebyshev) = A
 
+function ultraconv_nonequal(A, B)
+    a=order(A); b=order(B)
+    d=domain(A)
+    if -0.5 ≤ b-a ≤ 1 # -0.5, 0.5 or 1
+            return ConcreteConversion(A,B)
+    elseif b-a > 1
+        r = b:-1:a+1
+        v = [ConcreteConversion(Ultraspherical(i-1,d), Ultraspherical(i,d)) for i in r]
+        if !(last(r) ≈ a+1)
+            vlast = ConcreteConversion(A, Ultraspherical(last(r)-1, d))
+            v2 = Union{eltype(v), typeof(vlast)}[v; vlast]
+        else
+            v2 = v
+        end
+        bwsum = isapproxinteger(b-a) ? (0, 2length(v)) : (0,ℵ₀)
+        return ConversionWrapper(TimesOperator(v2, bwsum, (ℵ₀,ℵ₀), bwsum), A, B)
+    else
+        throw(ArgumentError("please implement $A → $B"))
+    end
+end
 
 function Conversion(A::Ultraspherical,B::Ultraspherical)
     @assert domain(A) == domain(B)
@@ -174,20 +194,7 @@ function Conversion(A::Ultraspherical,B::Ultraspherical)
     if b==a
         return ConversionWrapper(Operator(I,A))
     elseif isapproxinteger(b-a) || isapproxhalfoddinteger(b-a)
-        if -0.5 ≤ b-a ≤ 1 # -0.5, 0.5 or 1
-            return ConcreteConversion(A,B)
-        elseif b-a > 1
-            r = b:-1:a+1
-            v = [ConcreteConversion(Ultraspherical(i-1,d), Ultraspherical(i,d)) for i in r]
-            if !(last(r) ≈ a+1)
-                vlast = ConcreteConversion(A, Ultraspherical(last(r)-1, d))
-                v2 = Union{eltype(v), typeof(vlast)}[v; vlast]
-            else
-                v2 = v
-            end
-            bwsum = isapproxinteger(b-a) ? (0, 2length(v)) : (0,ℵ₀)
-            return ConversionWrapper(TimesOperator(v2, bwsum, (ℵ₀,ℵ₀), bwsum), A, B)
-        end
+        return ultraconv_nonequal(A, B)
     end
     throw(ArgumentError("please implement $A → $B"))
 end

test/JacobiTest.jl

@@ -75,25 +75,25 @@ include("testutils.jl")
                     @testset for J in (J1, NormalizedPolynomialSpace(J1))
                         g = Fun(f, J)
                         if !any(isnan, coefficients(g))
-                            @test Conversion(C, J) * f ≈ g
+                            @test @inferred(Conversion(C, J) * f) ≈ g
                         end
                     end
                 end
                 # Jacobi(-0.5, -0.5) to NormalizedJacobi(-0.5, -0.5) can have a NaN in the (1,1) index
                 # if the normalization is not carefully implemented,
                 # so this test checks that this is not the case
-                g = Conversion(Jacobi(-0.5,-0.5), NormalizedJacobi(-0.5,-0.5)) * Fun(x->x^3, Jacobi(-0.5, -0.5))
+                g = @inferred Conversion(Jacobi(-0.5,-0.5), NormalizedJacobi(-0.5,-0.5)) * Fun(x->x^3, Jacobi(-0.5, -0.5))
                 @test coefficients(g) ≈ coefficients(Fun(x->x^3, NormalizedChebyshev()))
             end
-            @testset "legendre" begin
+            @testset "Legendre" begin
                 f = Fun(x->x^2, Legendre(d))
                 C = space(f)
                 for J1 in (Jacobi(-0.5, -0.5, d), Legendre(d),
                                 Jacobi(0.5, 0.5, d), Jacobi(1, 1, d))
                     @testset for J in (J1, NormalizedPolynomialSpace(J1))
                         g = Fun(f, J)
                         if !any(isnan, coefficients(g))
-                            @test Conversion(C, J) * f ≈ g
+                            @test @inferred(Conversion(C, J) * f) ≈ g
                         end
                     end
                 end
@@ -106,7 +106,7 @@ include("testutils.jl")
                     @testset for J in (J1, NormalizedPolynomialSpace(J1))
                         g = Fun(f, J)
                         if !any(isnan, coefficients(g))
-                            @test Conversion(U, J) * f ≈ g
+                            @test @inferred(Conversion(U, J) * f) ≈ g
                         end
                     end
                 end
@@ -119,7 +119,7 @@ include("testutils.jl")
                     @testset for J in (J1,)
                         g = Fun(f, J)
                         if !any(isnan, coefficients(g))
-                            @test Conversion(U, J) * f ≈ g
+                            @test @inferred(Conversion(U, J) * f) ≈ g
                         end
                     end
                 end
@@ -129,14 +129,32 @@ include("testutils.jl")
         @testset for (S1, S2) in ((Legendre(), Jacobi(-0.5, -0.5)), (Jacobi(-0.5, -0.5), Legendre()),
                 (Legendre(), Jacobi(1.5, 1.5)), (Legendre(), Ultraspherical(0.5)),
                 (Ultraspherical(0.5), Legendre()))
-            @test Conversion(S1, S2) * Fun(x->x^4, S1) ≈ Fun(x->x^4, S2)
+            @test @inferred(Conversion(S1, S2) * Fun(x->x^4, S1)) ≈ Fun(x->x^4, S2)
         end
 
         C = Conversion(Jacobi(Chebyshev()), Ultraspherical(1))
-        @test C * Fun(x->x^2, Jacobi(-0.5, -0.5)) ≈ Fun(x->x^2, Ultraspherical(1))
+        @test @inferred(C * Fun(x->x^2, Jacobi(-0.5, -0.5))) ≈ Fun(x->x^2, Ultraspherical(1))
 
         C = Conversion(Ultraspherical(0.5), Jacobi(Chebyshev()))
-        @test C * Fun(x->x^2, Ultraspherical(0.5)) ≈ Fun(x->x^2, Jacobi(Chebyshev()))
+        @test @inferred(C * Fun(x->x^2, Ultraspherical(0.5))) ≈ Fun(x->x^2, Jacobi(Chebyshev()))
+
+        for Lb in -0.5:0.5:2, La in -0.5:0.5:2,
+            Mb in Lb:1:4, Ma in La:1:4
+                f = Fun(x -> x^4, Jacobi(Lb, La))
+                g = Fun(x -> x^4, Jacobi(Mb, Ma))
+                h = @inferred Conversion(Jacobi(Lb, La), Jacobi(Mb, Ma)) * f
+                @test space(h) == space(g)
+                @test h ≈ g
+        end
+
+        for Lb in -0.5:0.5:2, Mb in Lb:0.5:4
+            La = Lb; Ma = Mb
+            f = Fun(x -> x^4, Jacobi(Lb, La))
+            g = Fun(x -> x^4, Jacobi(Mb, Ma))
+            h = @inferred Conversion(Jacobi(Lb, La), Jacobi(Mb, Ma)) * f
+            @test space(h) == space(g)
+            @test h ≈ g
+        end
 
         @testset "inference tests" begin
             #= Note all cases are inferred as of now,