Improve type-inference in recursive Jacobi multiplication (#89)

* Improve type-inference in recursive Jacobi multiplication

* half integer tests

* Add inference test for HalfOddInt weights

Author: jishnub

Project.toml

@@ -1,10 +1,11 @@
 name = "ApproxFunSingularities"
 uuid = "f8fcb915-6b99-5be2-b79a-d6dbef8e6e7e"
-version = "0.3.16"
+version = "0.3.17"
 
 [deps]
 ApproxFunBase = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"
 ApproxFunOrthogonalPolynomials = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
+BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
@@ -25,6 +26,7 @@ ApproxFunBase = "0.8.56, 0.9.12"
 ApproxFunBaseTest = "0.1"
 ApproxFunOrthogonalPolynomials = "0.2.3, 0.3, 0.4, 0.5, 0.6"
 Aqua = "0.7"
+BlockBandedMatrices = "0.11, 0.12"
 DomainSets = "0.4, 0.5, 0.6"
 HalfIntegers = "1.5"
 IntervalSets = "0.5, 0.6, 0.7"

src/ApproxFunSingularities.jl

@@ -31,14 +31,16 @@ import ApproxFunBase: Fun, SumSpace, SubSpace, WeightSpace, NoSpace,
             dotu, components, promoterangespace, ∞,
             assert_integer, SpecialEvalPtType, isleftendpoint, isrightendpoint, evaluation_point,
             @calculus_operator, ConcreteConversion, InterlaceOperator_Diagonal, UnsetSpace,
-            choosedomainspace, mean
+            choosedomainspace, mean, bandwidthssum
 
 import ApproxFunOrthogonalPolynomials: order
 
 import IntervalSets: rightendpoint, leftendpoint, Domain
 
 import DomainSets: ChebyshevInterval, boundary, dimension
 
+using BlockBandedMatrices: blockbandwidths, subblockbandwidths
+
 import Base: convert, getindex, *, /, ^,
             show, sum, cumsum, complex, sqrt, abs, in, first, last,
             union, isapprox, zeros, one, length, ones, exp, log

src/JacobiWeight.jl

@@ -422,12 +422,22 @@ function Multiplication(f::Fun{<:JacobiWeight{<:ConstantSpace,<:IntervalOrSegmen
     elseif isapproxinteger(Sf.β) && Sf.β ≥ 1 && S.b >0
         # decrement β and multiply again
         M1 = ConcreteMultiplication(f.coefficients[1]*jacobiweight(1,0,d),S)
-        M1_out = Multiplication(jacobiweight(Sf.β-1,Sf.α,d), rangespace(M1)) * M1
+        M1_ = Multiplication(jacobiweight(Sf.β-1,Sf.α,d), rangespace(M1).space)
+        bw = bandwidths(M1) .+ bandwidths(M1_)
+        bbw = blockbandwidths(M1) .+ blockbandwidths(M1_)
+        sbbw = subblockbandwidths(M1) .+ subblockbandwidths(M1_)
+        ts = (size(M1, 1), size(M1_, 2))
+        M1_out = TimesOperator(Operator{eltype(M1)}[M1_, M1], bw, ts, bbw, sbbw)
         MultiplicationWrapper(f, M1_out, S)
     elseif isapproxinteger(Sf.α) && Sf.α ≥ 1 && S.a >0
         # decrement α and multiply again
         M2 = ConcreteMultiplication(f.coefficients[1]*jacobiweight(0,1,d),S)
-        M2_out = Multiplication(jacobiweight(Sf.β,Sf.α-1,d), rangespace(M2)) * M2
+        M2_ = Multiplication(jacobiweight(Sf.β,Sf.α-1,d), rangespace(M2).space)
+        bw = bandwidths(M2) .+ bandwidths(M2_)
+        bbw = blockbandwidths(M2) .+ blockbandwidths(M2_)
+        sbbw = subblockbandwidths(M2) .+ subblockbandwidths(M2_)
+        ts = (size(M2, 1), size(M2_, 2))
+        M2_out = TimesOperator(Operator{eltype(M2)}[M2_, M2], bw, ts, bbw, sbbw)
         MultiplicationWrapper(f, M2_out, S)
     else
         # default JacobiWeight
@@ -444,15 +454,18 @@ Multiplication(f::Fun{<:JacobiWeight{<:ConstantSpace,<:IntervalOrSegmentDomain}}
 function rangespace(M::ConcreteMultiplication{<:
                         JacobiWeight{<:ConstantSpace,<:IntervalOrSegmentDomain},<:Jacobi})
     S=domainspace(M)
-    if space(M.f).β==1
+    Sf = space(M.f)
+    zeroT = zero(Sf.β)
+    J = if Sf.β == 1
         # multiply by (1+x)
         Jacobi(S.b-1,S.a,domain(S))
-    elseif space(M.f).α == 1
+    elseif Sf.α == 1
         # multiply by (1-x)
         Jacobi(S.b,S.a-1,domain(S))
     else
         error("Not implemented")
     end
+    JacobiWeight(zeroT, zeroT, J)
 end
 
 bandwidths(::ConcreteMultiplication{<:

test/runtests.jl

@@ -3,7 +3,8 @@ module AFSTests
 using ApproxFunBase
 using ApproxFunBase: HeavisideSpace, PointSpace, ArraySpace, DiracSpace, PiecewiseSegment,
                         UnionDomain, resizedata!, CachedOperator, RaggedMatrix,
-                        Block, ∞, BandedBlockBandedMatrix, NoSpace
+                        Block, ∞, BandedBlockBandedMatrix, NoSpace, ConcreteMultiplication,
+                        MultiplicationWrapper
 using ApproxFunBaseTest: testbandedoperator, testtransforms, testfunctional,
                         testbandedblockbandedoperator
 using ApproxFunOrthogonalPolynomials
@@ -109,9 +110,57 @@ end
         xg = Fun(x->x*√(1-x^2), JacobiWeight(0.5, 0.5, Jacobi(1,1)))
         @test Multiplication(g) * Fun(NormalizedLegendre()) ≈ xg
 
-        f = Fun(x -> (1-x)^2, JacobiWeight(3,4,ConstantSpace(ChebyshevInterval())));
-        S = @inferred (f -> domainspace(Multiplication(f, Jacobi(1,1))))(f)
-        @test S == Jacobi(1,1)
+        genf(β,α) = Fun(x -> (1+x)^β * (1-x)^α, JacobiWeight(β,α,ConstantSpace(ChebyshevInterval())));
+
+        f = genf(0,2)
+        S = Jacobi(5,5)
+        d = domain(S)
+        # Multiplication(f, S) is inferred as a small Union
+        # We enumerate the possible types
+        T1 = typeof(Multiplication(genf(1,0), S)::ConcreteMultiplication)
+        T2 = typeof(Multiplication(f, S)::MultiplicationWrapper)
+        T3 = typeof(Multiplication(genf(1,0), Legendre())::MultiplicationWrapper)
+        @inferred Union{T1,T2,T3} Multiplication(f, S)
+
+        dsp(f, S) = domainspace(Multiplication(f, S))
+        rsp(f, S) = rangespace(Multiplication(f, S))
+        ds = if VERSION >= v"1.9"
+            @inferred dsp(f, S)
+        else
+            dsp(f, S)
+        end
+        @test ds == S
+        @test rsp(f, S) == Jacobi(5,3)
+
+        f = genf(4,3)
+        S = Jacobi(2,1)
+        @test rsp(f, S) == JacobiWeight(2,2,Legendre())
+
+        f = genf(half(Odd(3)), half(Odd(3)))
+        S = Jacobi(2,0)
+        @test @inferred(((f,S) -> domainspace(@inferred Multiplication(f, S)))(f,S)) == S
+        @test Multiplication(f, S) * Fun(S) ≈ Fun(x->x*(1-x^2)^(3/2), JacobiWeight(3/2, 3/2, S))
+
+        S = Jacobi(half(Odd(1)), half(Odd(3)))
+        @test @inferred(domainspace(@inferred Multiplication(f,S))) == S
+
+        @testset for β in 0:5, α in 0:5
+            f = genf(β, α)
+            g = Fun(x->(1+x)^2 * (1+x)^β * (1-x)^α)
+            @testset for b in 0:8, a in 0:8
+                S = Jacobi(b,a)
+                w = Fun(x->(1+x)^2, S)
+                M = Multiplication(f, S)
+                @test domainspace(M) == S
+                if b >= β && a >= α
+                    @test rangespace(M) == Jacobi(b-β, a-α)
+                end
+                if b < β && a < α
+                    @test rangespace(M) == JacobiWeight(β-b, α-a, Legendre())
+                end
+                @test Multiplication(f) * w ≈ M * w ≈ g
+            end
+        end
     end
 
     @testset "Derivative" begin