Use polynomial transforms from FastTransforms (#226)

* use poly transforms from FastTransforms

* add link to issue in test

* coefficient conversions to normalized polynomials

* short-circuit for same space

* disable ultra transforms

* reuse chebyshev-integer ultraspherical methods

* ultra2ultra for non-integer cases

* don't delete extra line in Chebyshev

* don't delete extra lines in Ultraspherical

* reduce infnorm tolerance in jacobi test

Author: jishnub

Project.toml

@@ -28,7 +28,7 @@ BlockArrays = "0.14, 0.15, 0.16"
 BlockBandedMatrices = "0.10, 0.11"
 DomainSets = "0.5, 0.6"
 FastGaussQuadrature = "0.4, 0.5"
-FastTransforms = "0.12, 0.13, 0.14, 0.15"
+FastTransforms = "0.12, 0.13, 0.14"
 FillArrays = "0.11, 0.12, 0.13"
 IntervalSets = "0.5, 0.6, 0.7"
 LazyArrays = "0.22"

src/Spaces/Jacobi/jacobitransform.jl

@@ -36,28 +36,3 @@ icjt(P, cfs, ::Val{false}) = P * cfs
 function *(P::JacobiITransformPlan, cfs::AbstractVector)
     P.ichebplan * icjt(P.icjtplan, cfs, Val(inplace(P)))
 end
-
-
-function coefficients(f::AbstractVector{T}, a::Jacobi, b::Chebyshev) where T
-    if domainscompatible(a, b) && !(isapproxinteger_addhalf(a.a) && isapproxinteger_addhalf(a.b))
-        jac2cheb(f, strictconvert(T,a.a), strictconvert(T,a.b))
-    else
-        defaultcoefficients(f,a,b)
-    end
-end
-function coefficients(f::AbstractVector{T}, a::Chebyshev, b::Jacobi) where T
-    isempty(f) && return f
-    if domainscompatible(a, b) && !(isapproxinteger_addhalf(b.a) && isapproxinteger_addhalf(b.b))
-        cheb2jac(f, strictconvert(T,b.a), strictconvert(T,b.b))
-    else
-        defaultcoefficients(f,a,b)
-    end
-end
-
-function coefficients(f::AbstractVector,a::Jacobi,b::Jacobi)
-    if domainscompatible(a, b) && !(isapproxinteger(a.a-b.a) && isapproxinteger(a.b-b.b))
-        jac2jac(f,a.a,a.b,b.a,b.b)
-    else
-        defaultcoefficients(f,a,b)
-    end
-end

src/Spaces/Spaces.jl

@@ -45,3 +45,131 @@ function getindex(C::ConcreteConversion{<:HeavisideSpace,
         k::Integer,j::Integer)
     k ≤ dimension(domainspace(C)) && j==k ? one(eltype(C)) : zero(eltype(C))
 end
+
+# Fast conversion between common bases using FastTransforms
+
+# The chebyshev-ultraspherical transforms are currently very slow,
+# see https://github.com/JuliaApproximation/FastTransforms.jl/issues/204
+# We therefore go through hoops to only call these for non-integral Ultraspherical orders
+
+function _changepolybasis(v::StridedVector{T},
+        C::MaybeNormalized{<:Chebyshev{<:ChebyshevInterval}},
+        U::MaybeNormalized{<:Ultraspherical{<:Any,<:ChebyshevInterval}},
+        ) where {T<:AbstractFloat}
+
+    normcheb = C isa NormalizedPolynomialSpace
+    normultra = U isa NormalizedPolynomialSpace
+    Uc = normultra ? _stripnorm(U) : U
+    Cc = normcheb ? _stripnorm(C) : C
+    if order(U) == 1
+        v = normcheb ? ApproxFunBase.mul_coefficients(ConcreteConversion(C, Cc), v) : v
+        vc = ultraconversion(v)
+        normultra ? ApproxFunBase.mul_coefficients!(ConcreteConversion(Uc, U), vc) : vc
+    elseif order(U) isa Union{Integer, StaticInt}
+        coefficients(v, C, Ultraspherical(1,domain(U)), U)
+    else
+        cheb2ultra(v, strictconvert(T, order(U)); normcheb, normultra)
+    end
+end
+function _changepolybasis(v::StridedVector{T},
+        U::MaybeNormalized{<:Ultraspherical{<:Any,<:ChebyshevInterval}},
+        C::MaybeNormalized{<:Chebyshev{<:ChebyshevInterval}}) where {T<:AbstractFloat}
+
+    normultra = U isa NormalizedPolynomialSpace
+    normcheb = C isa NormalizedPolynomialSpace
+    Uc = normultra ? _stripnorm(U) : U
+    Cc = normcheb ? _stripnorm(C) : C
+    if order(U) == 1
+        v = normultra ? ApproxFunBase.mul_coefficients(ConcreteConversion(U, Uc), v) : v
+        vc = ultraiconversion(v)
+        normcheb ? ApproxFunBase.mul_coefficients!(ConcreteConversion(Cc, C), vc) : vc
+    elseif order(U) isa Union{Integer, StaticInt}
+        coefficients(v, U, Ultraspherical(1,domain(U)), C)
+    else
+        ultra2cheb(v, strictconvert(T, order(U)); normultra, normcheb)
+    end
+end
+function _changepolybasis(v::StridedVector{T},
+        U1::MaybeNormalized{<:Ultraspherical{<:Any,<:ChebyshevInterval}},
+        U2::MaybeNormalized{<:Ultraspherical{<:Any,<:ChebyshevInterval}},
+        ) where {T<:AbstractFloat}
+
+    if isinteger(order(U1) - order(U2))
+        defaultcoefficients(v, U1, U2)
+    else
+        norm1 = U1 isa NormalizedPolynomialSpace
+        norm2 = U2 isa NormalizedPolynomialSpace
+        ultra2ultra(v, strictconvert(T, order(U1)), strictconvert(T, order(U2)); norm1, norm2)
+    end
+end
+
+function _changepolybasis(v::StridedVector{T},
+        C::MaybeNormalized{<:Chebyshev{<:ChebyshevInterval}},
+        J::MaybeNormalized{<:Jacobi{<:ChebyshevInterval}},
+        ) where {T<:AbstractFloat}
+
+    normcheb = C isa NormalizedPolynomialSpace
+    normjac = J isa NormalizedPolynomialSpace
+    Jc = _stripnorm(J)
+    if Jc.a == 0 && Jc.b == 0
+        cheb2leg(v; normcheb, normleg = normjac)
+    else
+        cheb2jac(v, strictconvert(T,Jc.a), strictconvert(T,Jc.b); normcheb, normjac)
+    end
+end
+function _changepolybasis(v::StridedVector{T},
+        J::MaybeNormalized{<:Jacobi{<:ChebyshevInterval}},
+        C::MaybeNormalized{<:Chebyshev{<:ChebyshevInterval}},
+        ) where {T<:AbstractFloat}
+
+    normcheb = C isa NormalizedPolynomialSpace
+    normjac = J isa NormalizedPolynomialSpace
+    Jc = _stripnorm(J)
+    if Jc.a == 0 && Jc.b == 0
+        leg2cheb(v; normcheb, normleg = normjac)
+    else
+        jac2cheb(v, strictconvert(T,Jc.a), strictconvert(T,Jc.b); normcheb, normjac)
+    end
+end
+function _changepolybasis(v::StridedVector{T},
+        U::MaybeNormalized{<:Ultraspherical{<:Any,<:ChebyshevInterval}},
+        J::MaybeNormalized{<:Jacobi{<:ChebyshevInterval}},
+        ) where {T<:AbstractFloat}
+
+    normultra = U isa NormalizedPolynomialSpace
+    normjac = J isa NormalizedPolynomialSpace
+    Jc = _stripnorm(J)
+    ultra2jac(v, strictconvert(T,order(U)), strictconvert(T,Jc.a), strictconvert(T,Jc.b);
+        normultra, normjac)
+end
+function _changepolybasis(v::StridedVector{T},
+        J::MaybeNormalized{<:Jacobi{<:ChebyshevInterval}},
+        U::MaybeNormalized{<:Ultraspherical{<:Any,<:ChebyshevInterval}},
+        ) where {T<:AbstractFloat}
+
+    normjac = J isa NormalizedPolynomialSpace
+    normultra = U isa NormalizedPolynomialSpace
+    Jc = _stripnorm(J)
+    jac2ultra(v, strictconvert(T,Jc.a), strictconvert(T,Jc.b), strictconvert(T,order(U));
+        normultra, normjac)
+end
+function _changepolybasis(v::StridedVector{T},
+        J1::MaybeNormalized{<:Jacobi{<:ChebyshevInterval}},
+        J2::MaybeNormalized{<:Jacobi{<:ChebyshevInterval}},
+        ) where {T<:AbstractFloat}
+
+    norm1 = J1 isa NormalizedPolynomialSpace
+    norm2 = J2 isa NormalizedPolynomialSpace
+    J1c = _stripnorm(J1)
+    J2c = _stripnorm(J2)
+    jac2jac(v, strictconvert(T,J1c.a), strictconvert(T,J1c.b), strictconvert(T,J2c.a), strictconvert(T,J2c.b);
+        norm1, norm2)
+end
+_changepolybasis(v, a, b) = defaultcoefficients(v, a, b)
+
+function coefficients(f::AbstractVector{T},
+        a::MaybeNormalized{<:Union{Chebyshev,Ultraspherical,Jacobi}},
+        b::MaybeNormalized{<:Union{Chebyshev,Ultraspherical,Jacobi}}) where T
+    spacescompatible(a,b) && return f
+    _changepolybasis(f, a, b)
+end

src/Spaces/Ultraspherical/UltrasphericalOperators.jl

@@ -190,8 +190,10 @@ function Conversion(A::Ultraspherical,B::Ultraspherical)
     throw(ArgumentError("please implement $A → $B"))
 end
 
-maxspace_rule(A::Ultraspherical,B::Ultraspherical) = order(A) > order(B) ? A : B
-
+function maxspace_rule(A::Ultraspherical, B::Ultraspherical)
+    isapproxinteger(order(A) - order(B)) || return NoSpace()
+    order(A) > order(B) ? A : B
+end
 
 function getindex(M::ConcreteConversion{<:Chebyshev,U,T},
         k::Integer,j::Integer) where {T, U<:Ultraspherical{<:Union{Integer, StaticInt}}}
@@ -292,27 +294,6 @@ function _conversion_rule(a::Ultraspherical, b::Ultraspherical)
     end
 end
 
-
-function coefficients(g::AbstractVector,sp::Ultraspherical{<:Union{Integer,StaticInt}},C::Chebyshev)
-    domainscompatible(C,sp) || throw(ArgumentError("domains don't match: $(domain(sp)) and $(domain(C))"))
-    if order(sp) == 1
-        ultraiconversion(g)
-    else
-        # do one at a time
-        coefficients(g,sp,Ultraspherical(1,domain(sp)),C)
-    end
-end
-function coefficients(g::AbstractVector,C::Chebyshev,sp::Ultraspherical)
-    domainscompatible(C,sp) || throw(ArgumentError("domains don't match: $(domain(C)) and $(domain(sp))"))
-    if order(sp) == 1
-        ultraconversion(g)
-    else
-        # do one at a time
-        coefficients(g,C,Ultraspherical(1,domain(sp)),sp)
-    end
-end
-
-
 # TODO: include in getindex to speed up
 function Integral(sp::Chebyshev{DD}, m::Number) where {DD<:IntervalOrSegment}
     assert_integer(m)

test/ChebyshevTest.jl

@@ -378,7 +378,7 @@ using ApproxFunOrthogonalPolynomials: forwardrecurrence
                     sp in Any[_sp, NormalizedPolynomialSpace(_sp)]
                 d = domain(sp)
                 f = Fun(sp, c)
-                for ep in [leftendpoint, rightendpoint],
+                @testset for ep in [leftendpoint, rightendpoint],
                         ev in [ApproxFunBase.ConcreteEvaluation, Evaluation]
                     E = @inferred ev(sp, ep, 0)
                     @test E[2:4] ≈ E[1:4][2:end]

test/JacobiTest.jl

@@ -344,9 +344,12 @@ using Static
         f = Fun(x->cospi(1000x))
         g = Fun(f,Legendre())
         h = Fun(g,Chebyshev())
-        @test norm(f.coefficients-h.coefficients,Inf) < 1000eps()
-        @time h = Fun(h,Legendre())
-        @test norm(g.coefficients-h.coefficients,Inf) < 10000eps()
+        @test norm(coefficients(f) - coefficients(h), Inf) < 1000eps()
+        @time j = Fun(h,Legendre())
+        # The accuracy in the following may be improved,
+        # See https://github.com/JuliaApproximation/FastTransforms.jl/issues/201
+        # On FastTransforms v0.15, an upper bound of 60000eps() works
+        @test norm(coefficients(g) - coefficients(j), Inf) < 10000eps()
     end
 
     @testset "conversion for non-compatible paramters" begin