Add tests for Polynomial spaces

Project.toml

@@ -39,22 +39,32 @@ DSP = "0.6, 0.7"
 DomainSets = "0.5, 0.6"
 DualNumbers = "0.6.2"
 FFTW = "0.3, 1"
+FastGaussQuadrature = "0.4, 0.5"
+FastTransforms = "0.12, 0.13, 0.14, 0.15.1"
 FillArrays = "0.11, 0.12, 0.13, 1"
+HalfIntegers = "1.5"
 InfiniteArrays = "0.11, 0.12"
 InfiniteLinearAlgebra = "0.5, 0.6"
 Infinities = "0.1"
 IntervalSets = "0.5, 0.6, 0.7"
 LazyArrays = "0.20, 0.21, 0.22, 1"
 LowRankApprox = "0.2, 0.3, 0.4, 0.5"
+OddEvenIntegers = "0.1.8"
 SimpleTraits = "0.9"
 SpecialFunctions = "0.10, 1.0, 2"
+Static = "0.8"
 StaticArrays = "0.12, 1.0"
 julia = "1.6"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c"
+HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721"
 Infinities = "e1ba4f0e-776d-440f-acd9-e1d2e9742647"
+OddEvenIntegers = "8d37c425-f37a-4ca2-9b9d-a61bc06559d2"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 
 [targets]
-test = ["Aqua", "Random", "Infinities"]
+test = ["Aqua", "Random", "FastTransforms", "FastGaussQuadrature", "HalfIntegers", "Infinities", "OddEvenIntegers", "Static"]

test/AFOP/Chebyshev/Chebyshev.jl

@@ -0,0 +1,258 @@
+
+export Chebyshev, NormalizedChebyshev
+
+"""
+`Chebyshev()` is the space spanned by the Chebyshev polynomials
+```
+    T_0(x),T_1(x),T_2(x),…
+```
+where `T_k(x) = cos(k*acos(x))`.  This is the default space
+as there exists a fast transform and general smooth functions on `[-1,1]`
+can be easily resolved.
+"""
+struct Chebyshev{D<:Domain,R} <: PolynomialSpace{D,R}
+    domain::D
+    function Chebyshev{D,R}(d) where {D,R}
+        isempty(d) && throw(ArgumentError("Domain cannot be empty"))
+        new(d)
+    end
+    Chebyshev{D,R}() where {D,R} = new(strictconvert(D, ChebyshevInterval()))
+end
+
+Chebyshev(d::Domain) = Chebyshev{typeof(d),real(prectype(d))}(d)
+Chebyshev() = Chebyshev(ChebyshevInterval())
+Chebyshev(d) = Chebyshev(Domain(d))
+const NormalizedChebyshev{D<:Domain,R} = NormalizedPolynomialSpace{Chebyshev{D, R}, D, R}
+NormalizedChebyshev() = NormalizedPolynomialSpace(Chebyshev())
+NormalizedChebyshev(d) = NormalizedPolynomialSpace(Chebyshev(d))
+
+function Base.getproperty(S::Chebyshev{<:Any,<:Any},v::Symbol)
+    if v==:b || v==:a
+        -0.5
+    else
+        getfield(S,v)
+    end
+end
+
+normalization(::Type{T}, sp::Chebyshev, k::Int) where T = T(π)/(2-FastTransforms.δ(k,0))
+
+Space(d::SegmentDomain) = Chebyshev(d)
+_norm(x) = norm(x)
+_norm(x::Real) = abs(x) # this preserves integers, whereas norm returns a float
+function Space(d::AbstractInterval)
+    a,b = endpoints(d)
+    d2 = if isinf(_norm(a)) && isinf(_norm(b))
+        Line(d)
+    elseif isinf(_norm(a)) || isinf(_norm(b))
+        Ray(d)
+    else
+        d
+    end
+    Chebyshev(d2)
+end
+
+
+setdomain(S::Chebyshev, d::Domain) = Chebyshev(d)
+
+ones(::Type{T}, S::Chebyshev) where {T<:Number} = Fun(S,fill(one(T),1))
+ones(S::Chebyshev) = Fun(S,fill(1.0,1))
+
+function first(f::Fun{<:Chebyshev})
+    n = ncoefficients(f)
+    n == 0 && return zero(cfstype(f))
+    n == 1 && return f.coefficients[1]
+    foldr(-,coefficients(f))
+end
+
+last(f::Fun{<:Chebyshev}) = reduce(+,coefficients(f))
+
+spacescompatible(a::Chebyshev,b::Chebyshev) = domainscompatible(a,b)
+hasfasttransform(::Chebyshev) = true
+supportsinplacetransform(::Chebyshev{<:Domain{T},T}) where {T<:AbstractFloat} = true
+
+
+## Transform
+
+transform(::Chebyshev,vals::AbstractVector,plan) = plan*vals
+itransform(::Chebyshev,cfs::AbstractVector,plan) = plan*cfs
+plan_transform(::Chebyshev,vals::AbstractVector) = plan_chebyshevtransform(vals)
+plan_itransform(::Chebyshev,cfs::AbstractVector) = plan_ichebyshevtransform(cfs)
+plan_transform!(::Chebyshev, vals::AbstractVector) = plan_chebyshevtransform!(vals)
+plan_itransform!(::Chebyshev, cfs::AbstractVector) = plan_ichebyshevtransform!(cfs)
+
+## Evaluation
+
+clenshaw(sp::Chebyshev, c::AbstractVector, x::AbstractArray) =
+    clenshaw(c,x,ClenshawPlan(promote_type(eltype(c),eltype(x)),sp,length(c),length(x)))
+
+clenshaw(::Chebyshev,c::AbstractVector,x) = chebyshev_clenshaw(c,x)
+
+clenshaw(c::AbstractVector,x::AbstractVector,plan::ClenshawPlan{S,V}) where {S<:Chebyshev,V}=
+    clenshaw(c,collect(x),plan)
+
+#TODO: This modifies x, which is not threadsafe
+clenshaw(c::AbstractVector,x::Vector,plan::ClenshawPlan{<:Chebyshev}) = chebyshev_clenshaw(c,x,plan)
+
+function clenshaw(c::AbstractMatrix{T},x::T,plan::ClenshawPlan{S,T}) where {S<:Chebyshev,T<:Number}
+    bk=plan.bk
+    bk1=plan.bk1
+    bk2=plan.bk2
+
+    m,n=size(c) # m is # of coefficients, n is # of funs
+
+    @inbounds for i = 1:n
+        bk1[i] = zero(T)
+        bk2[i] = zero(T)
+    end
+    x = 2x
+
+    @inbounds for k=m:-1:2
+        for j=1:n
+            ck = c[k,j]
+            bk[j] = muladd(x,bk1[j],ck - bk2[j])
+        end
+        bk2, bk1, bk = bk1, bk, bk2
+    end
+
+    x = x/2
+    @inbounds for i = 1:n
+        ce = c[1,i]
+        bk[i] = muladd(x,bk1[i],ce - bk2[i])
+    end
+
+    bk
+end
+
+function clenshaw(c::AbstractMatrix{T},x::AbstractVector{T},plan::ClenshawPlan{S,T}) where {S<:Chebyshev,T<:Number}
+    bk=plan.bk
+    bk1=plan.bk1
+    bk2=plan.bk2
+
+    m,n=size(c) # m is # of coefficients, n is # of funs
+
+    @inbounds for i = 1:n
+        x[i] = 2x[i]
+        bk1[i] = zero(T)
+        bk2[i] = zero(T)
+    end
+
+    @inbounds for k=m:-1:2
+        for j=1:n
+            ck = c[k,j]
+            bk[j] = muladd(x[j],bk1[j],ck - bk2[j])
+        end
+        bk2, bk1, bk = bk1, bk, bk2
+    end
+
+
+    @inbounds for i = 1:n
+        x[i] = x[i]/2
+        ce = c[1,i]
+        bk[i] = muladd(x[i],bk1[i],ce - bk2[i])
+    end
+
+    bk
+end
+
+# overwrite x
+
+function clenshaw!(c::AbstractVector,x::AbstractVector,plan::ClenshawPlan{S,V}) where {S<:Chebyshev,V}
+    N,n = length(c),length(x)
+
+    if isempty(c)
+        fill!(x,0)
+        return x
+    end
+
+    bk=plan.bk
+    bk1=plan.bk1
+    bk2=plan.bk2
+
+    @inbounds for i = 1:n
+        x[i] = 2x[i]
+        bk1[i] = zero(V)
+        bk2[i] = zero(V)
+    end
+
+    @inbounds for k = N:-1:2
+        ck = c[k]
+        for i = 1:n
+            bk[i] = muladd(x[i],bk1[i],ck - bk2[i])
+        end
+        bk2, bk1, bk = bk1, bk, bk2
+    end
+
+    ce = c[1]
+    @inbounds for i = 1:n
+        x[i] = x[i]/2
+        x[i] = muladd(x[i],bk1[i],ce-bk2[i])
+    end
+
+    x
+end
+
+## Calculus
+
+
+# diff T -> U, then convert U -> T
+integrate(f::Fun{Chebyshev{D,R}}) where {D<:IntervalOrSegment,R} =
+    Fun(f.space,fromcanonicalD(f,0)*ultraint!(ultraconversion(f.coefficients)))
+differentiate(f::Fun{Chebyshev{D,R}}) where {D<:IntervalOrSegment,R} =
+    Fun(f.space,1/fromcanonicalD(f,0)*ultraiconversion(ultradiff(f.coefficients)))
+
+
+
+
+
+## Multivariate
+
+# determine correct parameter to have at least
+# N point
+_padua_length(N) = Int(cld(-3+sqrt(1+8N),2))
+
+function squarepoints(::Type{T}, N) where T
+    pts = paduapoints(T, _padua_length(N))
+    n = size(pts,1)
+    ret = Array{SVector{2,T}}(undef, n)
+    @inbounds for k=1:n
+        ret[k] = SVector{2,T}(pts[k,1],pts[k,2])
+    end
+    ret
+end
+
+points(S::TensorSpace{<:Tuple{<:Chebyshev{<:ChebyshevInterval},<:Chebyshev{<:ChebyshevInterval}}}, N) =
+    squarepoints(real(prectype(S)), N)
+
+function points(S::TensorSpace{<:Tuple{<:Chebyshev,<:Chebyshev}},N)
+    T = real(prectype(S))
+    pts = squarepoints(T, N)
+    pts .= fromcanonical.(Ref(domain(S)), pts)
+    pts
+end
+
+plan_transform(S::TensorSpace{<:Tuple{<:Chebyshev,<:Chebyshev}},v::AbstractVector) =
+    plan_paduatransform!(v,Val{false})
+
+transform(S::TensorSpace{<:Tuple{<:Chebyshev,<:Chebyshev}},v::AbstractVector,
+        plan=plan_transform(S,v)) = plan*copy(v)
+
+plan_itransform(S::TensorSpace{<:Tuple{<:Chebyshev,<:Chebyshev}},v::AbstractVector) =
+     plan_ipaduatransform!(eltype(v),sum(1:Int(block(S,length(v)))),Val{false})
+
+itransform(S::TensorSpace{<:Tuple{<:Chebyshev,<:Chebyshev}},v::AbstractVector,
+         plan=plan_itransform(S,v)) = plan*pad(v,sum(1:Int(block(S,length(v)))))
+
+
+#TODO: adaptive
+for op in (:(Base.sin),:(Base.cos))
+    @eval ($op)(f::ProductFun{<:Chebyshev,<:Chebyshev}) =
+        ProductFun(chebyshevtransform($op.(values(f))),space(f))
+end
+
+
+
+reverseorientation(f::Fun{<:Chebyshev}) =
+    Fun(Chebyshev(reverseorientation(domain(f))),alternatesign!(copy(f.coefficients)))
+
+
+include("ChebyshevOperators.jl")