updates for new FastTransforms

.travis.yml

@@ -12,6 +12,8 @@ julia:
 matrix:
   allow_failures:
     - julia: nightly
+    - julia: "1.3"
+
 notifications:
   email: false
 codecov: true

Project.toml

@@ -28,11 +28,11 @@ BlockBandedMatrices = "0.5, 0.6"
 DomainSets = "0.1"
 FFTW = "0.3, 1"
 FastGaussQuadrature = "0.3.2, 0.4"
-FastTransforms = "0.5, 0.6"
+FastTransforms = "0.6.1"
 FillArrays = "0.7, 0.8"
 IntervalSets = "0.3.1"
 Reexport = "0.2"
-SpecialFunctions = "0.7, 0.8"
+SpecialFunctions = "0.7, 0.8, 0.9"
 julia = "1"
 
 [extras]

appveyor.yml

@@ -15,6 +15,7 @@ platform:
 matrix:
   allow_failures:
   - julia_version: nightly
+  - julia_version: 1.3
   - platform: x86
 
 branches:

src/ApproxFunOrthogonalPolynomials.jl

@@ -74,7 +74,7 @@ import LinearAlgebra: BlasInt, BlasFloat, norm, ldiv!, mul!, det, eigvals, dot,
 
 import FastTransforms: ChebyshevTransformPlan, IChebyshevTransformPlan, plan_chebyshevtransform,
                         plan_chebyshevtransform!, plan_ichebyshevtransform, plan_ichebyshevtransform!,
-                        pochhammer, lgamma
+                        pochhammer, lgamma, chebyshevtransform!, ichebyshevtransform!
 
 import BlockBandedMatrices: blockbandwidths, subblockbandwidths
 
@@ -94,51 +94,11 @@ import SpecialFunctions: sinpi, cospi, airy, besselh,
                     atan, acos, asin, erfc, inv
 
 points(d::IntervalOrSegmentDomain{T},n::Integer; kind::Int=1) where {T} =
-     fromcanonical.(Ref(d), chebyshevpoints(float(real(eltype(T))), n; kind=kind))  # eltype to handle point
-bary(v::AbstractVector{Float64}, d::IntervalOrSegmentDomain, x::Float64) = bary(v,tocanonical(d,x))
+    fromcanonical.(Ref(d), chebyshevpoints(float(real(eltype(T))), n; kind=kind))  # eltype to handle point
+bary(v::AbstractVector{Float64},d::IntervalOrSegmentDomain,x::Float64) = bary(v,tocanonical(d,x))
+
+include("bary.jl")
 
-function FastTransforms.chebyshevtransform!(X::AbstractMatrix{T};kind::Integer=1) where T<:fftwNumber		
-    if kind == 1		
-        if size(X) == (1,1)		
-            X		
-        else		
-            X=r2r!(X,REDFT10)		
-            X[:,1]/=2;X[1,:]/=2;		
-            lmul!(1/(size(X,1)*size(X,2)),X)		
-        end		
-    elseif kind == 2		
-        if size(X) == (1,1)		
-            X		
-        else		
-            X=r2r!(X,REDFT00)		
-            lmul!(1/((size(X,1)-1)*(size(X,2)-1)),X)		
-            X[:,1]/=2;X[:,end]/=2		
-            X[1,:]/=2;X[end,:]/=2		
-            X		
-        end		
-    end		
- end		
-
-function FastTransforms.ichebyshevtransform!(X::AbstractMatrix{T};kind::Integer=1) where T<:fftwNumber		
-    if kind == 1		
-        if size(X) == (1,1)		
-            X		
-        else		
-            X[1,:]*=2;X[:,1]*=2		
-            X = r2r(X,REDFT01)		
-            lmul!(1/4,X)		
-        end		
-    elseif kind == 2		
-        if size(X) == (1,1)		
-            X		
-        else		
-            X[1,:]*=2;X[end,:]*=2;X[:,1]*=2;X[:,end]*=2		
-            X=chebyshevtransform!(X;kind=kind)		
-            X[1,:]*=2;X[end,:]*=2;X[:,1]*=2;X[:,end]*=2		
-            lmul!((size(X,1)-1)*(size(X,2)-1)/4,X)		
-        end		
-    end		
-end
 
 include("ultraspherical.jl")
 include("Domains/Domains.jl")

src/Spaces/Jacobi/jacobitransform.jl

@@ -10,7 +10,7 @@ JacobiTransformPlan(chebplan::CPLAN, cjtplan::CJT) where {CPLAN,CJT} =
     JacobiTransformPlan{eltype(chebplan),CPLAN,CJT}(chebplan, cjtplan)
 
 plan_transform(S::Jacobi, v::AbstractVector) =
-    JacobiTransformPlan(plan_transform(Chebyshev(), v), plan_icjt(v, S.a, S.b))
+    JacobiTransformPlan(plan_transform(Chebyshev(), v), plan_cheb2jac(v, S.a, S.b))
 *(P::JacobiTransformPlan, vals::AbstractVector) = P.cjtplan*(P.chebplan*vals)
 
 
@@ -25,29 +25,29 @@ JacobiITransformPlan(chebplan::CPLAN, cjtplan::CJT) where {CPLAN,CJT} =
 
 
 plan_itransform(S::Jacobi, v::AbstractVector) =
-    JacobiITransformPlan(plan_itransform(Chebyshev(), v), plan_cjt(v, S.a, S.b))
+    JacobiITransformPlan(plan_itransform(Chebyshev(), v), plan_jac2cheb(v, S.a, S.b))
 *(P::JacobiITransformPlan, cfs::AbstractVector) = P.ichebplan*(P.icjtplan*cfs)
 
 
-function coefficients(f::AbstractVector,a::Jacobi,b::Chebyshev)
+function coefficients(f::AbstractVector{T}, a::Jacobi, b::Chebyshev) where T
     if domain(a) == domain(b) && (!isapproxinteger(a.a-0.5) || !isapproxinteger(a.b-0.5))
-        cjt(f,a.a,a.b)
+        jac2cheb(f, convert(T,a.a), convert(T,a.b))
     else
         defaultcoefficients(f,a,b)
     end
 end
-function coefficients(f::AbstractVector,a::Chebyshev,b::Jacobi)
+function coefficients(f::AbstractVector{T}, a::Chebyshev, b::Jacobi) where T
     isempty(f) && return f
     if domain(a) == domain(b) && (!isapproxinteger(b.a-0.5) || !isapproxinteger(b.b-0.5))
-        icjt(f,b.a,b.b)
+        cheb2jac(f, convert(T,b.a), convert(T,b.b))
     else
         defaultcoefficients(f,a,b)
     end
 end
 
 function coefficients(f::AbstractVector,a::Jacobi,b::Jacobi)
     if domain(a) == domain(b) && (!isapproxinteger(a.a-b.a) || !isapproxinteger(a.b-b.b))
-        jjt(f,a.a,a.b,b.a,b.b)
+        jac2jac(f,a.a,a.b,b.a,b.b)
     else
         defaultcoefficients(f,a,b)
     end

src/Spaces/Ultraspherical/Ultraspherical.jl

@@ -58,7 +58,7 @@ end
 function UltrasphericalPlan(λ::Number,vals)
     if λ == 0.5
         cp = plan_transform(Chebyshev(),vals)
-        c2lp = FastTransforms.th_cheb2legplan(eltype(vals),length(vals))
+        c2lp = plan_cheb2leg(eltype(vals),length(vals))
         UltrasphericalPlan{typeof(cp),typeof(c2lp)}(cp,c2lp)
     else
         error("Not implemented")
@@ -68,7 +68,7 @@ end
 function UltrasphericalIPlan(λ::Number,cfs)
     if λ == 0.5
         cp=plan_itransform(Chebyshev(),cfs)
-        c2lp=FastTransforms.th_leg2chebplan(eltype(cfs),length(cfs))
+        c2lp=plan_leg2cheb(eltype(cfs),length(cfs))
         UltrasphericalIPlan{typeof(cp),typeof(c2lp)}(cp,c2lp)
     else
         error("Not implemented")

src/bary.jl

@@ -0,0 +1,54 @@
+## Barycentric formula
+
+export bary,barysum
+
+function bary(v::AbstractVector{Float64},pts::AbstractVector{Float64},x::Float64)
+  n=length(v)
+  @assert n == length(pts)
+
+
+  retd = 1/(2*(x-pts[1]))
+  retn = v[1]*retd
+
+  for i = 2:2:n-1
+    @inbounds cd = 1/(x-pts[i])
+    retd -= cd
+    @inbounds retn -= v[i]*cd
+  end
+
+  for i = 3:2:n-1
+    @inbounds cd = 1/(x-pts[i])
+    retd += cd
+    @inbounds retn += v[i]*cd
+  end
+
+  cd = .5*(-1.)^(n-1)/(x-pts[n])
+  retd += cd
+  retn += v[n]*cd
+
+  retn/retd
+end
+
+
+bary(v::AbstractVector{Float64},x::Float64)=bary(v,chebyshevpoints(length(v);kind=2),x)
+
+function randomadaptivebary(f::Function)
+  r=rand()
+  fr=f(r)
+  err=1.
+  logn=1
+
+  tol=200eps()
+
+  vals=Float64[]
+
+  while err > tol
+      n=2^logn + 1
+      pts=chebyshevpoints(n;kind=2)
+      vals=f(pts)
+      err=abs(bary(vals,pts,r)-fr)/n
+      logn+=1
+  end
+
+  vals
+end

test/JacobiTest.jl

@@ -1,8 +1,8 @@
 using ApproxFunOrthogonalPolynomials, ApproxFunBase, Test, SpecialFunctions, LinearAlgebra
-    import ApproxFunBase: testbandedbelowoperator, testbandedoperator, testspace, testtransforms, Vec,
-                        maxspace, NoSpace, hasconversion, testfunctional,
-                        reverseorientation, ReverseOrientation
-    import ApproxFunOrthogonalPolynomials: jacobip
+import ApproxFunBase: testbandedbelowoperator, testbandedoperator, testspace, testtransforms, Vec,
+                    maxspace, NoSpace, hasconversion, testfunctional,
+                    reverseorientation, ReverseOrientation
+import ApproxFunOrthogonalPolynomials: jacobip
 
 @testset "Jacobi" begin
     @testset "Basic" begin

test/MultivariateTest.jl

@@ -2,7 +2,7 @@ using ApproxFunBase, ApproxFunOrthogonalPolynomials, LinearAlgebra, SpecialFunct
 import ApproxFunBase: testbandedblockbandedoperator, testraggedbelowoperator, factor, Block, cfstype,
                     blocklengths, block, tensorizer, Vec, ArraySpace, ∞,
                     testblockbandedoperator
-import ApproxFunOrthogonalPolynomials: chebyshevtransform
+import ApproxFunOrthogonalPolynomials: chebyshevtransform                    
 
 @testset "Multivariate" begin
     @testset "Square" begin