Avoid recursion in Chebyshev -> Ultraspherical (#191)

* Avoid recursion in Chebyshev() -> Ultraspherical()

* precompute bandwidths

Author: jishnub

src/Spaces/Jacobi/JacobiOperators.jl

@@ -178,7 +178,7 @@ function Conversion(L::Jacobi,M::Jacobi)
         end
     end
 
-    error("Implement for $L → $M")
+    throw(ArgumentError("please implement $A → $B"))
 end
 
 bandwidths(::ConcreteConversion{<:Jacobi,<:Jacobi}) = (0,1)

src/Spaces/Ultraspherical/Ultraspherical.jl

@@ -36,6 +36,11 @@ setdomain(S::Ultraspherical,d::Domain) = Ultraspherical(order(S),d)
 
 
 convert(::Type{Ultraspherical{T,D,R}}, S::Ultraspherical{T,D,R}) where {T,D,R} = S
+convert(::Type{Ultraspherical{A,D,R}}, S::Ultraspherical{B,D,R}) where {A,B,D,R} =
+    Ultraspherical{A,D,R}(convert(A, order(S)), domain(S))
+
+promote_rule(::Type{Ultraspherical{A,D,R}}, ::Type{Ultraspherical{B,D,R}}) where {A,B,D,R} =
+    Ultraspherical{promote_type(A,B),D,R}
 
 
 canonicalspace(S::Ultraspherical) = Chebyshev(domain(S))

src/Spaces/Ultraspherical/UltrasphericalOperators.jl

@@ -124,28 +124,33 @@ end
 
 ## Conversion Operator
 
+isequalhalf(x) = x == 0.5
+isequalhalf(::Integer) = false
 
-function Conversion(A::Chebyshev,B::Ultraspherical)
-    if order(B) ≤ 1
-        ConcreteConversion(A,B)
-    else
-        d=domain(A)
-        US=Ultraspherical(order(B)-1,d)
-        ConversionWrapper(TimesOperator(
-            ConcreteConversion(US,B), Conversion(A,US)))
+function Conversion(A::Chebyshev, B::Ultraspherical)
+    mB = order(B)
+    d=domain(A)
+    dB = domain(B)
+    d == dB || throw(ArgumentError("domains must be identical"))
+    if isequalhalf(mB) || mB == 1
+        return ConcreteConversion(A,B)
+    elseif (isinteger(mB) || isapproxinteger_addhalf(mB)) && mB > 0
+        r = mB:-1:(isinteger(mB) ? 2 : 1)
+        v = [ConcreteConversion(Ultraspherical(i-1, d), Ultraspherical(i,d)) for i in r]
+        U = domainspace(last(v))
+        CAU = ConcreteConversion(A, U)
+        v2 = Union{eltype(v), typeof(CAU)}[v; CAU]
+        bwsum = (0, (isapproxinteger(mB) ? 2length(v2) : ℵ₀))
+        return ConversionWrapper(TimesOperator(v2, bwsum, (ℵ₀,ℵ₀), bwsum))
     end
+    throw(ArgumentError("please implement $A → $B"))
 end
 
-
-isequalhalf(x) = x == 0.5
-isequalhalf(::Integer) = false
-
 function Conversion(A::Ultraspherical,B::Chebyshev)
-    if isequalhalf(order(A))
-        ConcreteConversion(A,B)
-    else
-        error("Not implemented")
+    if isequalhalf(order(A)) && domain(A) == domain(B)
+        return ConcreteConversion(A,B)
     end
+    throw(ArgumentError("please implement $A → $B"))
 end
 
 
@@ -154,27 +159,26 @@ maxspace_rule(A::Ultraspherical,B::Chebyshev) = A
 
 function Conversion(A::Ultraspherical,B::Ultraspherical)
     a=order(A); b=order(B)
+    d=domain(A)
+    dB = domain(B)
+    d == dB || throw(ArgumentError("domains must be identical"))
     if b==a
-        ConversionWrapper(Operator(I,A))
+        return ConversionWrapper(Operator(I,A))
     elseif isapproxinteger(b-a) || isapproxinteger_addhalf(b-a)
         if -1 ≤ b-a ≤ 1 && (a,b) ≠ (2,1)
-            ConcreteConversion(A,B)
+            return ConcreteConversion(A,B)
         elseif b-a > 1
-            d=domain(A)
-            US=Ultraspherical(b-static(1),d)
             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))
                 v = [v; vlast]
             end
-            ConversionWrapper(TimesOperator(v))
-        else
-            throw(ArgumentError("Cannot convert from $A to $B"))
+            bwsum = (0, (isapproxinteger(b-a) ? 2length(v) : ℵ₀))
+            return ConversionWrapper(TimesOperator(v, bwsum, (ℵ₀,ℵ₀), bwsum))
         end
-    else
-        throw(ArgumentError("Cannot convert from $A to $B"))
     end
+    throw(ArgumentError("please implement $A → $B"))
 end
 
 maxspace_rule(A::Ultraspherical,B::Ultraspherical) = order(A) > order(B) ? A : B

test/UltrasphericalTest.jl

@@ -6,6 +6,9 @@ using LinearAlgebra
 using Static
 
 @verbose @testset "Ultraspherical" begin
+    @testset "promotion" begin
+        @inferred (() -> [Ultraspherical(1), Ultraspherical(2.0)])()
+    end
     @testset "identity fun" begin
         for d in (ChebyshevInterval(), 3..4, Segment(2, 5), Segment(1, 4im)), order in (1, 2, 0.5)
             s = Ultraspherical(order, d)
@@ -25,27 +28,27 @@ using Static
 
         # Conversions from Chebyshev to Ultraspherical should lead to a small union of types
         Tallowed = Union{
-            ApproxFunBase.ConcreteConversion{
-                Chebyshev{ChebyshevInterval{Float64}, Float64},
-                Ultraspherical{Int64, ChebyshevInterval{Float64}, Float64}, Float64},
-            ApproxFunBase.ConversionWrapper{TimesOperator{Float64, Tuple{Int64, Int64}}, Float64}};
+            typeof(Conversion(Chebyshev(), Ultraspherical(1))),
+            typeof(Conversion(Chebyshev(), Ultraspherical(2)))};
         @inferred Tallowed Conversion(Chebyshev(), Ultraspherical(1));
         @inferred Tallowed Conversion(Chebyshev(), Ultraspherical(2));
+        @inferred Conversion(Chebyshev(), Ultraspherical(static(2)));
+        @inferred Conversion(Chebyshev(), Ultraspherical(static(0.5)));
+        @inferred (() -> Conversion(Chebyshev(), Ultraspherical(2)))();
+        @inferred (() -> Conversion(Chebyshev(), Ultraspherical(0.5)))();
+        @inferred (() -> Conversion(Chebyshev(), Ultraspherical(2.5)))();
+
         # Conversions between Ultraspherical should lead to a small union of types
         Tallowed = Union{
-            ApproxFunBase.ConcreteConversion{
-                Ultraspherical{Int64, ChebyshevInterval{Float64}, Float64},
-                Ultraspherical{Int64, ChebyshevInterval{Float64}, Float64}, Float64},
-            ApproxFunBase.ConversionWrapper{
-                ConstantOperator{Float64,
-                    Ultraspherical{Int64, ChebyshevInterval{Float64}, Float64}}, Float64},
-                    ApproxFunBase.ConversionWrapper{TimesOperator{Float64, Tuple{Int64, Int64}}, Float64}};
+            typeof(Conversion(Ultraspherical(1), Ultraspherical(2))),
+            typeof(Conversion(Ultraspherical(1), Ultraspherical(3)))}
         @inferred Tallowed Conversion(Ultraspherical(1), Ultraspherical(2));
+        @inferred Tallowed Conversion(Ultraspherical(1), Ultraspherical(3));
 
         @inferred Conversion(Ultraspherical(static(1)), Ultraspherical(static(4)))
-        # these cases should be handled by constant-propagation
-        @inferred (() -> Conversion(Ultraspherical(static(1)), Ultraspherical(4)))()
-        @inferred (() -> Conversion(Ultraspherical(1), Ultraspherical(static(4))))()
+        @inferred (() -> Conversion(Ultraspherical(1), Ultraspherical(4)))()
+        @inferred (() -> Conversion(Ultraspherical(1.0), Ultraspherical(4.0)))();
+        @inferred (() -> Conversion(Ultraspherical(1.0), Ultraspherical(3.5)))();
 
         for n in (2,5)
             C1 = Conversion(Chebyshev(), Ultraspherical(n))
@@ -62,19 +65,24 @@ using Static
         g = CLC * f
         @test g ≈ Fun(x->x^2, Chebyshev())
 
-        f = Fun(x->x^2, Chebyshev()) # Legendre
-        CCL = Conversion(Chebyshev(), Ultraspherical(0.5))
-        @test !isdiag(CCL)
-        g = CCL * f
-        @test g ≈ Fun(x->x^2, Ultraspherical(0.5))
-
-        f = Fun(x->x^2, Ultraspherical(0.5)) # Legendre
-        for n in (2.5, 3)
-            CLU = Conversion(Ultraspherical(0.5), Ultraspherical(n))
-            @test !isdiag(CLU)
-            g = CLU * f
+        for n in (0.5, 1, 1.5, 2)
+            f = Fun(x->x^2, Chebyshev())
+            CCL = Conversion(Chebyshev(), Ultraspherical(n))
+            @test !isdiag(CCL)
+            g = CCL * f
             @test g ≈ Fun(x->x^2, Ultraspherical(n))
         end
+
+        ff = x->x^2 +2x^3 + 3x^4
+        for n1 in (0.5, 1)
+            f = Fun(ff, Ultraspherical(n1))
+            for n2 in (2.5, 3)
+                CLU = Conversion(Ultraspherical(n1), Ultraspherical(n2))
+                @test !isdiag(CLU)
+                g = CLU * f
+                @test g ≈ Fun(ff, Ultraspherical(n1))
+            end
+        end
     end
 
     @testset "Normalized space" begin
@@ -110,12 +118,12 @@ using Static
             @test transform(S, v) ≈ transform(J, v)
         end
         @testset for T in (Float32, Float64), ET in (T, complex(T))
-            v = Array{ET}(undef, 10)
-            v2 = similar(v)
-            M = Array{ET}(undef, 10, 10)
-            M2 = similar(M)
-            A = Array{ET}(undef, 10, 10, 10)
-            A2 = similar(A)
+            v = Array{ET}(undef, 10);
+            v2 = similar(v);
+            M = Array{ET}(undef, 10, 10);
+            M2 = similar(M);
+            A = Array{ET}(undef, 10, 10, 10);
+            A2 = similar(A);
             @testset for d in ((), (0..1,)), order in (0.5, 0.7, 1.5, 1, 3)
                 U = Ultraspherical(order, d...)
                 NU = NormalizedPolynomialSpace(U)