type inference improvements in CosSpace/Fourier integrate (#83)

* type inference improvements in CosSpace/Fourier integrate

* Test for empty case

* Test against AFOP for CosSpace integrate

Author: jishnub

Project.toml

@@ -19,6 +19,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 AbstractFFTs = "0.5, 1"
 ApproxFunBase = "0.8.6"
 ApproxFunBaseTest = "0.1"
+ApproxFunOrthogonalPolynomials = "0.6"
 Aqua = "0.5"
 BandedMatrices = "0.16, 0.17"
 DomainSets = "0.3, 0.4, 0.5, 0.6"
@@ -33,9 +34,10 @@ julia = "1.6"
 
 [extras]
 ApproxFunBaseTest = "a931bfaf-0cfd-4a5c-b69c-bf2eed002b43"
+ApproxFunOrthogonalPolynomials = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ApproxFunBaseTest", "Aqua", "Test", "SpecialFunctions"]
+test = ["ApproxFunBaseTest", "ApproxFunOrthogonalPolynomials", "Aqua", "Test", "SpecialFunctions"]

src/calculus.jl

@@ -68,20 +68,25 @@ linesum(f::Fun{CosSpace{DD,RR}}) where {DD<:PeriodicSegment,RR} =
 
 
 function integrate(f::Fun{CS}) where CS<:CosSpace
-    if isa(domain(f),Circle)
-        error("Integrate not implemented for CosSpace on Circle")
-    else  # Probably periodic itnerval, drop constant term if zero
+    d = domain(f)
+    tinf = if d isa Circle
+        throw(ArgumentError("Integrate not implemented for CosSpace on Circle"))
+    else  # Probably periodic interval, drop constant term if zero
         tol=1E-14 #TODO: smart tolerance.  Here relative is a bit tricky
                   # since this is called by Fourier integrate
-        if abs(f.coefficients[1])<tol
+        if isempty(f.coefficients) || abs(f.coefficients[1])<tol
             integrate(Fun(f,space(f)|(2:∞)))
-        else
-            d=domain(f)
-            @assert isa(d,PeriodicSegment)
+        elseif d isa PeriodicSegment
             x=Fun(identity, Interval(d))
-            (f.coefficients[1]*x)⊕integrate(Fun(f,space(f)|(2:∞)))
+            fconstant = (f.coefficients[1]*x)
+            fperiodic = integrate(Fun(f,space(f)|(2:∞)))
+            SS = SumSpace(space(fconstant), space(fperiodic))
+            Fun(SS, ApproxFunBase.interlace(coefficients(fconstant), coefficients(fperiodic)))
+        else
+            throw(ArgumentError("not implemented"))
         end
     end
+    chop(tinf)
 end
 
 function integrate(f::Fun{SS}) where SS<:SinSpace
@@ -98,9 +103,13 @@ for OP in (:differentiate,:integrate)
     @eval $OP(f::Fun{Fourier{D,R},T}) where {T,D<:Circle,R} = $OP(Fun(f,Laurent))
 end
 
-integrate(f::Fun{Fourier{D,R},T}) where {T,D<:PeriodicSegment,R} =
-    integrate(component(f,2))⊕integrate(component(f,1))
-
+function integrate(f::Fun{Fourier{D,R},T}) where {T,D<:PeriodicSegment,R}
+    g1 = integrate(component(f,2))
+    g2 = integrate(component(f,1))
+    # g1 ⊕ g2, but we construct the space directly to improve type-stability
+    SS = SumSpace(space(g1), space(g2))
+    Fun(SS, ApproxFunBase.interlace(coefficients(g1), coefficients(g2)))
+end
 
 
 

test/runtests.jl

@@ -4,6 +4,7 @@ using ApproxFunBase: Block, EmptyDomain, UnionDomain
 using ApproxFunBaseTest: testspace, testtransforms, testmultiplication, testraggedbelowoperator,
                     testbandedoperator, testblockbandedoperator, testbandedblockbandedoperator,
                     testcalculus, testfunctional
+using ApproxFunOrthogonalPolynomials
 using LinearAlgebra
 using SpecialFunctions
 _factorial(n) = gamma(n+1)
@@ -299,6 +300,12 @@ end
                 Fun(θ->sin(sin(θ)),Fourier()),Fun(θ->cos(θ)+cos(3θ),CosSpace()))
         @test norm(integrate(f)'-f)<10eps()
     end
+    @test iszero(integrate(Fun(CosSpace(), Float64[])))
+    f = Fun(θ->1+cos(θ), CosSpace())
+    g = integrate(f)
+    θ = 0.4
+    @test g(θ) ≈ θ + sin(θ)
+    @test Fun(g', space(f)) ≈ f
 end