Use matrix of points for Disk

Author: dlfivefifty

src/Cone/Cone.jl

@@ -246,18 +246,18 @@ rectspace(::DuffyCone) = NormalizedJacobi(0,1,Segment(1,0))*ZernikeDisk()
 
 blocklengths(d::DuffyCone) = blocklengths(rectspace(d))
 
-# M = N(N+1)/2
-# N*M == N^2(N+1)/2 == n
-# N = 1/3 (-1 + 1/(-1 + 27n + 3sqrt(3)sqrt(n*(-2 + 27 n)))^(1/3) + (-1 + 27n + 3sqrt(3)sqrt(n*(-2 + 27n)))^(1/3))
+# M = N*(2N-1)
+# N*N*(2N-1) == N^2*(2N-1) == n
+# N =1/6*(1 + 1/(1 + 54n + 6sqrt(3)sqrt(n + 27n^2))^(1/3) + (1 + 54n + 6sqrt(3)sqrt(n + 27n^2))^(1/3))
 function pointsize(::Cone, n) 
-    N = round(Int,1/3*(-1 + 1/(-1 + 27n + 3sqrt(3)sqrt(n*(-2 + 27n)))^(1/3) + (-1 + 27n + 3sqrt(3)sqrt(n*(-2 + 27n)))^(1/3)),RoundUp)
-    N, N*(N+1) ÷ 2
+    N = round(Int, 1/6*(1 + 1/(1 + 54n + 6sqrt(3)sqrt(n + 27n^2))^(1/3) + (1 + 54n + 6sqrt(3)sqrt(n + 27n^2))^(1/3)), RoundUp)
+    N, N, 2N-1
 end
 
 
 function points(d::Cone,n)
     N,M = pointsize(d,n)
-    pts = Array{float(eltype(d))}(undef,0)
+    pts = Array{ApproxFunBase.float(eltype(d))}(undef,0)
     a,b = rectspace(DuffyCone()).spaces
     for y in points(b,M), x in points(a,N)
         push!(pts,conemap(Vec(x...,y...)))
@@ -434,7 +434,13 @@ function LegendreConeTransformPlan(S::LegendreCone, v::AbstractVector{T}) where
 end
 
 
-*(P::LegendreConeTransformPlan, v::AbstractVector) = _coefficients(P.triangleplan, P.duffyplan*v, DuffyCone(), LegendreCone())
+function *(P::LegendreConeTransformPlan, v::AbstractVector) 
+    N,M = P.duffyplan.plan[1][2],P.duffyplan.plan[2][2]
+    V=reshape(v,N,M) 
+    cfs = P.duffyplan*copy(V)
+    duffy2legendrecone!(P.triangleplan,cfs)
+    fromtensor(LegendreCone(), cfs)
+end
 
 plan_transform(K::LegendreCone, v::AbstractVector) = LegendreConeTransformPlan(K, v)
 

src/Disk/Disk.jl

@@ -21,14 +21,10 @@ end
 # N = (1 + sqrt(1 + 8n)/4)
 function points(d::ZernikeDisk, n)
     a,b = rectspace(ZernikeDisk()).spaces
-    pts=Array{float(eltype(domain(d)))}(undef,0)
-    N,M = pointsize(d, n)
-
-    for y in points(b,M),
-        x in points(a,N)
-        push!(pts,polar(Vec(x...,y...)))
-    end
-    pts
+    M,N = pointsize(d, n)
+    p_b = points(b,N)
+    p_a = points(a,M)
+    polar.(Vec.(p_a, p_b'))
 end
 
 function evaluate(cfs::AbstractVector, Z::ZernikeDisk, xy) 
@@ -183,24 +179,22 @@ struct ZernikeDiskTransformPlan{DUF,CHEB}
     disk2cxf::CHEB
 end
 
-function ZernikeDiskTransformPlan(S::ZernikeDisk, v::AbstractVector{T}) where T
-    n = length(v)
-    N = (1 + isqrt(1+8n)) ÷ 4
-    M = 2N-1
+function ZernikeDiskTransformPlan(S::ZernikeDisk, V::AbstractMatrix{T}) where T
+    N,M = size(V)
+    @assert M == 2N-1
     D = plan_transform!(rectspace(ChebyshevDisk()), Array{T}(undef,N,M))
     N = (M-1)÷4+1
     4(N-1)+1 == M || (N += 1) # round up
     ZernikeDiskTransformPlan(D, CDisk2CxfPlan(N))
 end
 
-function *(P::ZernikeDiskTransformPlan, v::AbstractVector) 
-    N,M = P.cxfplan.plan[1][2],P.cxfplan.plan[2][2] 
-    V = P.cxfplan*reshape(copy(v),N,M)
+function *(P::ZernikeDiskTransformPlan, V::AbstractMatrix) 
+    V = P.cxfplan*copy(V)
     C = checkerboard(V)
     fromtensor(ZernikeDisk(),  P.disk2cxf\C)
 end
 
-plan_transform(K::ZernikeDisk, v::AbstractVector) = ZernikeDiskTransformPlan(K, v)
+plan_transform(K::ZernikeDisk, v::AbstractMatrix) = ZernikeDiskTransformPlan(K, v)
 
 struct ZernikeDiskITransformPlan{DUF,CHEB}
     icxfplan::DUF

test/test_cone.jl

@@ -2,7 +2,7 @@ using ApproxFun, MultivariateOrthogonalPolynomials, StaticArrays, FillArrays, Bl
 import ApproxFunBase: checkpoints, plan_transform!, fromtensor
 import MultivariateOrthogonalPolynomials: rectspace, totensor, duffy2legendreconic!, legendre2duffyconic!, c_plan_rottriangle, plan_transform,
                         c_execute_tri_hi2lo, c_execute_tri_lo2hi, duffy2legendrecone_column_J!, duffy2legendrecone!, legendre2duffycone!,
-                        duffy2legendreconic!, ConicTensorizer
+                        duffy2legendreconic!, ConicTensorizer, pointsize
 
 
 @testset "Conic" begin
@@ -126,6 +126,11 @@ end
         @test checkpoints(rs) isa Vector{SVector{3,Float64}}
     end
 
+    @testset "ConeTensorizer" begin
+        ts = MultivariateOrthogonalPolynomials.ConeTensorizer()
+        @test totensor(ts,1:10) == [1 2 3 5 6 7; 4 8 9 0 0 0; 10 0 0 0 0 0]
+    end
+
     @testset "DuffyCone" begin
         p = points(DuffyCone(), 10)
         @test p isa Vector{SVector{3,Float64}}
@@ -205,27 +210,33 @@ end
     @testset "LegendreConePlan" begin
         for N = 1:10
             n = N*sum(1:N)
-            @test round(Int,1/3*(-1 + 1/(-1 + 27n + 3sqrt(3)sqrt(n*(-2 + 27n)))^(1/3) + (-1 + 27n + 3sqrt(3)sqrt(n*(-2 + 27n)))^(1/3)),RoundUp) ≈ N
+            @test round(Int, 1/6*(1 + 1/(1 + 54n + 6sqrt(3)sqrt(n + 27n^2))^(1/3) + (1 + 54n + 6sqrt(3)sqrt(n + 27n^2))^(1/3)), RoundUp)
         end
         p = points(LegendreCone(),10)
-        @test length(p) == 18 == 3*sum(1:3)
+        @test length(p) == 12 == 2*2*3
         v = fill(1.0,length(p))
         n = length(v)
-        N = round(Int,1/3*(-1 + 1/(-1 + 27n + 3sqrt(3)sqrt(n*(-2 + 27n)))^(1/3) + (-1 + 27n + 3sqrt(3)sqrt(n*(-2 + 27n)))^(1/3)),RoundUp)
-        M = N*(N+1) ÷ 2
+        M,N,O = pointsize(Cone(),n)
+        
+
         D = plan_transform!(rectspace(DuffyCone()), reshape(v,N,M))
         N,M = D.plan[1][2],D.plan[2][2]
         V=reshape(v,N,M)
-        @test D*copy(V) ≈ [[1/sqrt(2) zeros(1,5)]; zeros(2,6)]
+        @test D*V ≈ [[1/sqrt(2) zeros(1,5)]; zeros(2,6)]
         T = Float64
         P = c_plan_rottriangle(N, zero(T), zero(T), zero(T))
         duffy2legendrecone!(P,V)
-        @test V ≈ [[1 zeros(1,2)]; zeros(1,3)]
-        @test fromtensor(DuffyConic(), V) ≈ [1; Zeros(5)]
+        @test V ≈ [[1/sqrt(2) zeros(1,5)]; zeros(2,6)]
+        @test fromtensor(LegendreCone(), V) ≈ [1/sqrt(2); Zeros(55)]
 
         v = fill(1.0,length(p))
-        @test plan_transform(LegendreConic(),v)*v ≈ [1; Zeros(3)]
+        P = plan_transform(LegendreCone(),v)
+        @test P*v ≈ [1/sqrt(2); Zeros(55)]
         @test v == fill(1.0,length(p))
+
+        p = points(LegendreCone(),20)    
+        v = fill(1.0,length(p))
+        P = plan_transform(LegendreCone(),v)
 
 
         for (m,k,ℓ) in ((0,0,0), )

test/test_disk.jl

@@ -141,23 +141,20 @@ end
     @testset "transform derivation" begin
         p = points(ZernikeDisk(),10)
         @test length(p) == 6
-        v = fill(1.0,length(p))
-        n = length(v)
-        N = (1 + isqrt(1+8n)) ÷ 4
-        M = 2N-1
-        D = plan_transform!(rectspace(ChebyshevDisk()), reshape(v,N,M))
+        @test size(p) == (2,3)
+        V = fill(1.0,size(p))
+        D = plan_transform!(rectspace(ChebyshevDisk()), V)
         N,M = D.plan[1][2],D.plan[2][2]
-        V=reshape(v,N,M)
         @test D*V ≈ [[1 zeros(1,2)]; zeros(1,3)]
         C = checkerboard(V)
         disk2cxf = CDisk2CxfPlan(size(C,1))
         C = disk2cxf\C
         @test C ≈ [[1/sqrt(2) zeros(1,4)]; zeros(1,5)]
         @test fromtensor(ZernikeDisk(), C) ≈ [1/sqrt(2); Zeros(5)]
 
-        v = fill(1.0,length(p))
+        v = fill(1.0,size(p))
         @test plan_transform(ZernikeDisk(),v)*v ≈ [1/sqrt(2); Zeros(5)]
-        @test v == fill(1.0,length(p))
+        @test v == fill(1.0,size(p))
     end
 
     @testset "ChebyshevDisk-ZernikeDisk" begin