getindex for Cos/SinSpace ConcreteEvaluation (#89)

* getindex for Cos/SinSpace ConcreteEvaluation

* Number arguments in SinSpace

* ConcreteEvaluation for Fourier

* ConcreteEvaluation for Laurent

* rename _eval to _conceval

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunFourier"
 uuid = "59844689-9c9d-51bf-9583-5b794ec66d30"
-version = "0.3.19"
+version = "0.3.20"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
@@ -17,7 +17,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 AbstractFFTs = "0.5, 1"
-ApproxFunBase = "0.8.9"
+ApproxFunBase = "0.8.16"
 ApproxFunBaseTest = "0.1"
 ApproxFunOrthogonalPolynomials = "0.6"
 Aqua = "0.5"

src/ApproxFunFourier.jl

@@ -35,7 +35,8 @@ import ApproxFunBase: Fun, SumSpace, SubSpace, NoSpace, IntervalOrSegment,
             reverseeven!, negateeven!, cfstype, alternatesign!, extremal_args,
             hesseneigvals, chebyshev_clenshaw, roots, EmptyDomain,
             chebmult_getindex, components, affine_setdiff, complexroots,
-            assert_integer, companion_matrix, InterlaceOperator_Diagonal
+            assert_integer, companion_matrix, InterlaceOperator_Diagonal,
+            evaluation_point, SpecialEvalPtType
 
 import BandedMatrices: bandwidths
 
@@ -261,6 +262,31 @@ function evaluate(f::AbstractVector,S::CosSpace,t)
     end
 end
 
+function _conceval(C, ::CosSpace, order, x::SpecialEvalPtType, m)
+    if iseven(order)
+        one(eltype(C))
+    else
+        zero(eltype(C))
+    end
+end
+function _conceval(C, ::CosSpace, order, x, m)
+    y = tocanonical(domain(C), x)
+    t = if iseven(order)
+        cos(m * y)
+    else
+        sin(m * y)
+    end
+    convert(eltype(C), t)
+end
+function getindex(C::ConcreteEvaluation{<:CosSpace{<:PeriodicSegment}}, k::Integer)
+    m = k - 1
+    order = C.order
+    L = period(domain(C))
+    t = _conceval(C, domainspace(C), C.order, evaluation_point(C), m)
+    s = mod(order, 4) ∈ (1,2) ? -1 : 1
+    r = s * (m * 2pi/L)^order * t
+    convert(eltype(C), r)
+end
 
 points(sp::SinSpace,n)=points(domain(sp),2n+2)[2:n+1]
 
@@ -292,7 +318,31 @@ itransform(sp::SinSpace,vals::AbstractVector,plan) = plan*vals
 
 evaluate(f::AbstractVector,S::SinSpace,t) = sineshaw(f,tocanonical(S,t))
 
-
+function _conceval(C, ::SinSpace, order, x::SpecialEvalPtType, k)
+    if iseven(order)
+        zero(eltype(C))
+    else
+        one(eltype(C))
+    end
+end
+function _conceval(C, ::SinSpace, order, x, m)
+    y = tocanonical(domain(C), x)
+    t = if iseven(order)
+        sin(m * y)
+    else
+        cos(m * y)
+    end
+    convert(eltype(C), t)
+end
+function getindex(C::ConcreteEvaluation{<:SinSpace{<:PeriodicSegment}}, k::Integer)
+    m = k
+    order = C.order
+    L = period(domain(C))
+    t = _conceval(C, domainspace(C), C.order, evaluation_point(C), m)
+    s = mod(order, 4) ∈ (0,1) ? 1 : -1
+    r = s * (m * 2pi/L)^order * t
+    convert(eltype(C), r)
+end
 
 ## Laurent space
 """

src/FourierOperators.jl

@@ -379,3 +379,30 @@ for TYP in (:Fourier,:Laurent,:CosSpace,:SinSpace,:Taylor)
         end
     end
 end
+
+Evaluation(F::Fourier{<:PeriodicSegment}, x, order) = ConcreteEvaluation(F, x, order)
+function _conceval(C, ::Fourier, order, x::SpecialEvalPtType, m, k)
+    if (iseven(order) && isodd(k)) || (isodd(order) && iseven(k))
+        one(eltype(C))
+    else
+        zero(eltype(C))
+    end
+end
+function _conceval(C, ::Fourier, order, x, m, k)
+    y = tocanonical(domain(C), x)
+    t = if (iseven(order) && isodd(k)) || (isodd(order) && iseven(k))
+        cos(m * y)
+    else
+        sin(m * y)
+    end
+    convert(eltype(C), t)
+end
+function getindex(C::ConcreteEvaluation{<:Fourier{<:PeriodicSegment}}, k::Integer)
+    m = k ÷ 2
+    order = C.order
+    L = period(domain(C))
+    t = _conceval(C, domainspace(C), C.order, evaluation_point(C), m, k)
+    s = (-1)^((order + isodd(k))÷2)
+    r = s * (m * 2pi/L)^order * t
+    convert(eltype(C), r)
+end

src/LaurentOperators.jl

@@ -315,3 +315,24 @@ function Conversion(A::Laurent{DD,RR},B::Laurent{DD,RR}) where {DD,RR}
         InterlaceOperator(Diagonal([Matrix(I,1,1),PermutationOperator([2,1])]))
     ,A,B))
 end
+
+Evaluation(F::Laurent{<:PeriodicSegment}, x, order) = ConcreteEvaluation(F, x, order)
+function _conceval(C, ::Laurent, order, x::SpecialEvalPtType, m, k)
+    one(eltype(C))
+end
+flip_even(k) = isodd(k) ? 1 : -1
+function _conceval(C, ::Laurent, order, x, m, k)
+    y = tocanonical(domain(C), x)
+    s = flip_even(k)
+    t = cis(s * m * y)
+    convert(eltype(C), t)
+end
+function getindex(C::ConcreteEvaluation{<:Laurent{<:PeriodicSegment}}, k::Integer)
+    m = k ÷ 2
+    order = C.order
+    L = period(domain(C))
+    t = _conceval(C, domainspace(C), C.order, evaluation_point(C), m, k)
+    s = flip_even(k)
+    r = (im * s * m * 2pi/L)^order * t
+    convert(eltype(C), r)
+end

test/runtests.jl

@@ -110,6 +110,35 @@ end
 
     ## Bug in multiplicaiton
     @test Fun(SinSpace(),Float64[])^2 == Fun(SinSpace(),Float64[])
+
+    @testset "Evaluation" begin
+        @testset "CosSpace" begin
+            @testset for sp in [CosSpace(), CosSpace(1..2)]
+                f = Fun(sp, Float64[1:8;])
+                @test ldirichlet(sp) * f ≈ f(leftendpoint(domain(f)))
+                @test rdirichlet(sp) * f ≈ f(rightendpoint(domain(f)))
+                @test (lneumann(sp) * f)(leftendpoint(domain(f))) ≈ f'(leftendpoint(domain(f))) atol=1e-10
+                @test (rneumann(sp) * f)(rightendpoint(domain(f))) ≈ f'(rightendpoint(domain(f))) atol=1e-10
+                pt = 1.2
+                @test Evaluation(sp, pt, 0) * f ≈ f(pt)
+                @test Evaluation(sp, pt, 1) * f ≈ f'(pt)
+                @test Evaluation(sp, pt, 2) * f ≈ f''(pt)
+            end
+        end
+        @testset "SinSpace" begin
+            @testset for sp in [SinSpace(), SinSpace(1..2)]
+                f = Fun(sp, Float64[1:8;])
+                @test (ldirichlet(sp) * f)(leftendpoint(domain(f))) ≈ f(leftendpoint(domain(f))) atol=1e-10
+                @test (rdirichlet(sp) * f)(leftendpoint(domain(f))) ≈ f(rightendpoint(domain(f))) atol=1e-10
+                @test lneumann(sp) * f ≈ f'(leftendpoint(domain(f)))
+                @test rneumann(sp) * f ≈ f'(rightendpoint(domain(f)))
+                pt = 1.2
+                @test Evaluation(sp, pt, 0) * f ≈ f(pt)
+                @test Evaluation(sp, pt, 1) * f ≈ f'(pt)
+                @test Evaluation(sp, pt, 2) * f ≈ f''(pt)
+            end
+        end
+    end
 end
 
 
@@ -252,6 +281,20 @@ end
         @test coefficients([4; zeros(4); [1, 2]; zeros(4); 3], Fourier(0..6pi), Fourier(0..2pi)) ≈ [4,1,2,3]
         @test coefficients([4; zeros(6); [1, 2]; zeros(6); 3], Fourier(0..8pi), Fourier(0..2pi)) ≈ [4,1,2,3]
     end
+
+    @testset "Evaluation" begin
+        @testset for sp in [Fourier(), Fourier(1..2)]
+            f = Fun(sp, Float64[1:8;])
+            @test ldirichlet(sp) * f ≈ f(leftendpoint(domain(f)))
+            @test rdirichlet(sp) * f ≈ f(rightendpoint(domain(f)))
+            @test (lneumann(sp) * f)(leftendpoint(domain(f))) ≈ f'(leftendpoint(domain(f))) atol=1e-10
+            @test (rneumann(sp) * f)(rightendpoint(domain(f))) ≈ f'(rightendpoint(domain(f))) atol=1e-10
+            pt = 1.2
+            @test Evaluation(sp, pt, 0) * f ≈ f(pt)
+            @test Evaluation(sp, pt, 1) * f ≈ f'(pt)
+            @test Evaluation(sp, pt, 2) * f ≈ f''(pt)
+        end
+    end
 end
 
 @testset "Laurent" begin
@@ -271,6 +314,20 @@ end
     ## Diagonal Derivative
     D = @inferred Derivative(Laurent())
     @test isdiag(D)
+
+    @testset "Evaluation" begin
+        @testset for sp in [Laurent(), Laurent(1..2)]
+            f = Fun(sp, Float64[1:8;])
+            @test ldirichlet(sp) * f ≈ f(leftendpoint(domain(f)))
+            @test rdirichlet(sp) * f ≈ f(rightendpoint(domain(f)))
+            @test (lneumann(sp) * f)(leftendpoint(domain(f))) ≈ f'(leftendpoint(domain(f))) atol=1e-10
+            @test (rneumann(sp) * f)(rightendpoint(domain(f))) ≈ f'(rightendpoint(domain(f))) atol=1e-10
+            pt = 1.2
+            @test Evaluation(sp, pt, 0) * f ≈ f(pt)
+            @test Evaluation(sp, pt, 1) * f ≈ f'(pt)
+            @test Evaluation(sp, pt, 2) * f ≈ f''(pt)
+        end
+    end
 end
 
 @testset "Circle" begin