Update transform

Author: dlfivefifty

src/Disk/Disk.jl

@@ -7,11 +7,15 @@ end
 ZernikeDisk(d::Disk{V}) where V = ZernikeDisk{V,complex(eltype(V))}(d)
 ZernikeDisk() = ZernikeDisk(Disk())
 
+spacescompatible(::ZernikeDisk, ::ZernikeDisk) = true
+
 @containsconstants ZernikeDisk
 
 points(K::ZernikeDisk, n::Integer) =
     fromcanonical.(Ref(K), points(ChebyshevDisk(), n))
 
+evaluate(cfs::AbstractVector, Z::ZernikeDisk, xy) = 
+    Fun(Fun(Z, cfs), ChebyshevDisk())(xy)
 
 struct ChebyshevDisk{V,T} <: Space{Disk{V},T}
     domain::Disk{V}
@@ -84,6 +88,33 @@ fromrectcfs(S, v) = fromtensor(S, checkerboard(totensor(rectspace(S),v)))
 evaluate(cfs::AbstractVector, S::ChebyshevDisk, x) = evaluate(torectcfs(S,cfs), rectspace(S), ipolar(x))
 
 
+function _coefficients(disk2cxf, v̂::AbstractVector, ::ChebyshevDisk, ::ZernikeDisk)
+    F = totensor(ChebyshevDisk(), v̂)
+    F *= sqrt(convert(T,π))
+    n = size(F,1)
+    F̌ = disk2cxf \ pad(F,n,4n-3)
+    trivec(F̌)
+end
+
+function _coefficients(disk2cxf, v::AbstractVector, ::ZernikeDisk,  ::ChebyshevDisk)
+    F̌ = tridevec(v)
+    F̌ /= sqrt(convert(T,π))
+    n = size(F,1)
+    F = disk2cxf \ pad(F̌,n,4n-3)
+    fromtensor(ChebyshevDisk(), F)
+end
+
+function coefficients(cfs::AbstractVector, CD::ChebyshevDisk, ZD::ZernikeDisk)
+    c = totensor(ChebyshevDisk(), cfs)            # TODO: wasteful
+    n = size(c,1)
+    _coefficients(CDisk2CxfPlan(n), cfs, CD, ZD)
+end
+
+function coefficients(cfs::AbstractVector, ZD::ZernikeDisk, CD::ChebyshevDisk)
+    c = tridevec(cfs)            # TODO: wasteful
+    n = size(c,1)
+    _coefficients(CDisk2CxfPlan(n), cfs, CD, ZD)
+end
 
 
 struct FastZernikeDiskTransformPlan{DUF,CHEB}
@@ -94,15 +125,14 @@ end
 function FastZernikeDiskTransformPlan(S::ZernikeDisk, v::AbstractVector{T}) where T
     # n = floor(Integer,sqrt(2length(v)) + 1/2)
     # v = Array{T}(undef, sum(1:n))
-    FastZernikeDiskTransformPlan(plan_transform(ChebyshevDisk(),v), CDisk2CxfPlan(n))
+    cxfplan = plan_transform(ChebyshevDisk(),v)
+    c = totensor(ChebyshevDisk(), cxfplan*v)            # TODO: wasteful
+    n = size(c,1)
+    FastZernikeDiskTransformPlan(cxfplan, CDisk2CxfPlan(n))
 end
 
-function *(P::FastZernikeDiskTransformPlan, v)
-    v̂ = P.cxfplan*v
-    F = tridevec_trans(v̂)
-    F̌ = P.disk2cheb \ F
-    trivec(tridenormalize!(F̌))
-end
+*(P::FastZernikeDiskTransformPlan, v::AbstractVector) = 
+    _coefficients(P.disk2cxf, P.cxfplan*v, ChebyshevDisk(), ZernikeDisk())
 
 plan_transform(K::ZernikeDisk, v::AbstractVector) = FastZernikeDiskTransformPlan(K, v)
 

src/c_transforms.jl

@@ -66,7 +66,7 @@ end
 CDisk2CxfPlan(n::Int) = CDisk2CxfPlan(c_plan_disk2cxf(n), n)
 
 function *(C::CDisk2CxfPlan, A::Matrix{Float64})
-    (size(A,1) == C.n+1 && size(A,2) == 4C.n+1) || throw(ArgumentError(A))
+    (size(A,1) == C.n+1 && size(A,2) == 4C.n-3) || throw(ArgumentError(A))
     B = copy(A)
     c_disk2cxf(C.plan, B)
     B

test/test_disk.jl

@@ -1,19 +1,10 @@
 using Revise, ApproxFun, MultivariateOrthogonalPolynomials
 import MultivariateOrthogonalPolynomials: checkerboard, icheckerboard, CDisk2CxfPlan
 import ApproxFunBase: totensor
-import ApproxFunOrthogonalPolynomials: jacobip, chebyshevt
+import ApproxFunOrthogonalPolynomials: jacobip
 
 
 
-
-
-f = (x,y) -> x*y+cos(y-0.1)+sin(x)+1; 
-ff = Fun(f, ChebyshevDisk(), 1000)
-@test ff(0.1,0.2) ≈ f(0.1,0.2)
-
-ff = Fun(f, ChebyshevDisk())
-@test ff(0.1,0.2) ≈ f(0.1,0.2)
-
 function chebydiskeval(c::AbstractMatrix{T}, r, θ) where T
     ret = zero(T)
     for j = 1:2:size(c,2), k=1:size(c,1)
@@ -28,66 +19,83 @@ function chebydiskeval(c::AbstractMatrix{T}, r, θ) where T
 end
 
 
-for (m,ℓ) in ((0,1),(0,2),(0,3), (1,0), (1,1), (1,2), (2,0), (2,1),(3,0),(3,4))
-    p1 = (r,θ) -> cos(m*θ) * cos((2ℓ+isodd(m))*acos(r))
-    f = Fun((x,y) -> p1(sqrt(x^2+y^2), atan(y,x)), ChebyshevDisk())
-    c = totensor(f.space, chop(f.coefficients,1E-12))    
-    @test c[ℓ+1, 2m+1] ≈ 1
-    @test f(0.1cos(0.2),0.1sin(0.2)) ≈ chebydiskeval(c, 0.1, 0.2)
+@testset "ChebyshevDisk" begin
+    f = (x,y) -> x*y+cos(y-0.1)+sin(x)+1; 
+    ff = Fun(f, ChebyshevDisk(), 1000)
+    @test ff(0.1,0.2) ≈ f(0.1,0.2)
+
+    ff = Fun(f, ChebyshevDisk())
+    @test ff(0.1,0.2) ≈ f(0.1,0.2)
+
 
-    if m > 0
-        p2 = (r,θ) -> sin(m*θ) * cos((2ℓ+isodd(m))*acos(r))
-        f = Fun((x,y) -> p2(sqrt(x^2+y^2), atan(y,x)), ChebyshevDisk())
+    for (m,ℓ) in ((0,1),(0,2),(0,3), (1,0), (1,1), (1,2), (2,0), (2,1),(3,0),(3,4))
+        p1 = (r,θ) -> cos(m*θ) * cos((2ℓ+isodd(m))*acos(r))
+        f = Fun((x,y) -> p1(sqrt(x^2+y^2), atan(y,x)), ChebyshevDisk())
         c = totensor(f.space, chop(f.coefficients,1E-12))    
-        @test c[ℓ+1, 2m] ≈ 1
+        @test c[ℓ+1, 2m+1] ≈ 1
         @test f(0.1cos(0.2),0.1sin(0.2)) ≈ chebydiskeval(c, 0.1, 0.2)
+
+        if m > 0
+            p2 = (r,θ) -> sin(m*θ) * cos((2ℓ+isodd(m))*acos(r))
+            f = Fun((x,y) -> p2(sqrt(x^2+y^2), atan(y,x)), ChebyshevDisk())
+            c = totensor(f.space, chop(f.coefficients,1E-12))    
+            @test c[ℓ+1, 2m] ≈ 1
+            @test f(0.1cos(0.2),0.1sin(0.2)) ≈ chebydiskeval(c, 0.1, 0.2)
+        end
     end
 end
 
-f = Fun((x,y) -> 1, ChebyshevDisk()); n = 1;
+
+@testset "disk transform" begin
+    f = Fun((x,y) -> 1, ChebyshevDisk()); n = 1;
     c = totensor(f.space, f.coefficients)
     P = CDisk2CxfPlan(n)
     d = P \ c
     @test d ≈ [1/sqrt(2)]
-f = Fun((x,y) -> y, ChebyshevDisk()); n = 2;
+    f = Fun((x,y) -> y, ChebyshevDisk()); n = 2;
     c = totensor(f.space, f.coefficients)
     c = pad(c, n, 4n-3)
     P = CDisk2CxfPlan(n)
     d = P \ c
-    @test d ≈  [0.0  0.5  0.0  0.0  0.0; 0.0  0.0  0.0  0.0  0.0]
+    @test d ≈  [0.0  0.5  0.0 0.0 0.0 ; 0.0  0.0  0.0 0.0 0.0]
 
-f = Fun((x,y) -> x, ChebyshevDisk()); n = 2;
+    f = Fun((x,y) -> x, ChebyshevDisk()); n = 2;
     c = totensor(f.space, f.coefficients)
     c = pad(c, n, 4n-3)
     P = CDisk2CxfPlan(n)
     d = P \ c
-    @test d ≈  [0.0  0.0  0.5  0.0  0.0; 0.0  0.0  0.0  0.0  0.0]
-
-
-# m,ℓ = (0,2);
-# m,ℓ = (1,1);
-# m,ℓ = (2,2);
-
-
-m,ℓ = (2,4);
-p1 = (r,θ) -> sqrt(2ℓ+2)*r^m*jacobip((ℓ-m)÷2, m, 0, 2r^2-1)*cos(m*θ)/sqrt(2π)
-f = Fun((x,y) -> p1(sqrt(x^2+y^2), atan(y,x)), ChebyshevDisk());
-    @test f(0.1*cos(0.2),0.1*sin(0.2)) ≈ p1(0.1,0.2)
-    c = totensor(f.space, chop(f.coefficients,1E-12))
-    @test chebydiskeval(c, 0.1, 0.2) ≈ p1(0.1,0.2)
-    
-
-    n = size(c,1); c = sqrt(2π)*pad(c,n,4n-3)
-    P = CDisk2CxfPlan(n)
-    d = P \ c
-    @test d[(ℓ-m)÷2, 2m+1] ≈ 0
-    @test d[(ℓ-m)÷2+1, 2m+1] ≈ 1
-
-    
-
-
-
-    
-
-
+    @test d ≈  [0.0  0.0  0.5 0.0 0.0; 0.0  0.0  0.0 0.0 0.0]
+
+
+    @testset "explicit polynomial" begin
+        p1 = (r,θ) -> sqrt(10)*r^2*(4r^2-3)*cos(2θ)/sqrt(π)
+        f = Fun((x,y) -> p1(sqrt(x^2+y^2), atan(y,x)), ChebyshevDisk());
+        @test f(0.1*cos(0.2),0.1*sin(0.2)) ≈ p1(0.1,0.2)
+        c = totensor(f.space, chop(f.coefficients,1E-12))
+        @test chebydiskeval(c, 0.1, 0.2) ≈ p1(0.1,0.2)
+        n = size(c,1); c = sqrt(π)*pad(c,n,4n-3)
+        P = CDisk2CxfPlan(n)
+        d = P \ c
+        @test d[1,4] ≈ 0 atol=1E-13
+        @test d[2,5] ≈ 1
+    end
+    @testset "different m and ℓ" begin
+        for (m,ℓ) in ((0,2), (1,1), (2,2), (2,6), (3,3))
+            p1 = (r,θ) -> sqrt(2ℓ+2)*r^m*jacobip((ℓ-m)÷2, 0, m, 2r^2-1)*cos(m*θ)/sqrt(π)
+            f = Fun((x,y) -> p1(sqrt(x^2+y^2), atan(y,x)), ChebyshevDisk());
+            @test f(0.1*cos(0.2),0.1*sin(0.2)) ≈ p1(0.1,0.2)
+            c = totensor(f.space, chop(f.coefficients,1E-12))
+            @test chebydiskeval(c, 0.1, 0.2) ≈ p1(0.1,0.2)
+
+            n = size(c,1); c = sqrt(π)*pad(c,n,4n-3)
+            P = CDisk2CxfPlan(n)
+            d = P \ c
+            @test norm(d[:, 1:2m]) ≈ 0 atol = 1E-13
+            @test norm(d[1:(ℓ-m)÷2, 2m+1]) ≈ 0 atol = 1E-13
+            @test d[(ℓ-m)÷2+1, 2m+1] ≈ 1
+            @test norm(d[(ℓ-m)÷2+2:end, 2m+1]) ≈ 0 atol = 1E-13
+            @test norm(d[:, 2m+2:end]) ≈ 0 atol = 1E-13
+        end
+    end
+end