update disk2cxf, add rectdisk2cheb

Author: MikaelSlevinsky

Project.toml

@@ -25,7 +25,7 @@ BinaryProvider = "0.5"
 DSP = "0.6"
 FFTW = "1"
 FastGaussQuadrature = "0.4"
-FastTransforms_jll = "0.3.3"
+FastTransforms_jll = "0.4.0"
 FillArrays = "0.8, 0.9"
 Reexport = "0.2"
 SpecialFunctions = "0.8, 0.9, 0.10"

README.md

@@ -19,7 +19,7 @@ julia> using FastTransforms, LinearAlgebra
 
 ## Fast orthogonal polynomial transforms
 
-The 29 orthogonal polynomial transforms are listed in `FastTransforms.kind2string.(0:28)`. Univariate transforms may be planned with the standard normalization or with orthonormalization. For multivariate transforms, the standard normalization may be too severe for floating-point computations, so it is omitted. Here are two examples:
+The 33 orthogonal polynomial transforms are listed in `FastTransforms.kind2string.(0:32)`. Univariate transforms may be planned with the standard normalization or with orthonormalization. For multivariate transforms, the standard normalization may be too severe for floating-point computations, so it is omitted. Here are two examples:
 
 ### The Chebyshev--Legendre transform
 

examples/disk.jl

@@ -37,8 +37,9 @@ r = [sinpi((N-n-0.5)/(2N)) for n in 0:N-1]
 # On the mapped tensor product grid, our function samples are:
 F = [f(r*cospi(θ), r*sinpi(θ)) for r in r, θ in θ]
 
-# We precompute a Zernike--Chebyshev×Fourier plan:
-P = plan_disk2cxf(F)
+# We precompute a (generalized) Zernike--Chebyshev×Fourier plan:
+α, β = 0, 0
+P = plan_disk2cxf(F, α, β)
 
 # And an FFTW Chebyshev×Fourier analysis plan on the disk:
 PA = plan_disk_analysis(F)
@@ -47,10 +48,57 @@ PA = plan_disk_analysis(F)
 U = P\(PA*F)
 
 # The Zernike coefficients are useful for integration. The integral of $f(x,y)$
-# over the disk should be $\pi/2$ by harmonicity. The coefficient of $Z_0^0$
+# over the disk should be $\pi/2$ by harmonicity. The coefficient of $Z_{0,0}$
 # multiplied by `√π` is:
 U[1, 1]*sqrt(π)
 
 # Using an orthonormal basis, the integral of $[f(x,y)]^2$ over the disk is
 # approximately the square of the 2-norm of the coefficients:
-norm(U)^2
+norm(U)^2 - π/(2*sqrt(2))*log1p(sqrt(2))
+
+# But there's more! Next, we repeat the experiment using the Dunkl-Xu
+# orthonormal polynomials supported on the rectangularized disk.
+N = 2N
+M = N
+
+# We analyze the function on an $N\times M$ mapped tensor product $xy$-grid defined by:
+# ```math
+# \begin{aligned}
+# x_n & = \cos\left(\frac{2n+1}{2N}\pi\right) = \sin\left(\frac{N-2n-1}{2N}\pi\right),\quad {\rm for} \quad 0 \le n < N,\quad{\rm and}\\
+# z_m & = \cos\left(\frac{2m+1}{2M}\pi\right) = \sin\left(\frac{M-2m-1}{2M}\pi\right),\quad {\rm for} \quad 0 \le m < M,\\
+# y_{n,m} & = \sqrt{1-x_n^2}z_m.
+# \end{aligned}
+# ```
+# Slightly more accuracy can be expected by using an auxiliary array:
+# ```math
+#   w_n = \sin\left(\frac{2n+1}{2N}\pi\right),\quad {\rm for} \quad 0 \le n < N,
+# ```
+# so that $y_{n,m} = w_nz_m$.
+#
+# The x grid
+w = [sinpi((n+0.5)/N) for n in 0:N-1]
+x = [sinpi((N-2n-1)/(2N)) for n in 0:N-1]
+
+# The z grid
+z = [sinpi((M-2m-1)/(2M)) for m in 0:M-1]
+
+# On the mapped tensor product grid, our function samples are:
+F = [f(x[n], w[n]*z) for n in 1:N, z in z]
+
+# We precompute a Dunkl-Xu--Chebyshev plan:
+P = plan_rectdisk2cheb(F, β)
+
+# And an FFTW Chebyshev² analysis plan on the rectangularized disk:
+PA = plan_rectdisk_analysis(F)
+
+# Its Dunkl-Xu coefficients are:
+U = P\(PA*F)
+
+# The Dunkl-Xu coefficients are useful for integration. The integral of $f(x,y)$
+# over the disk should be $\pi/2$ by harmonicity. The coefficient of $P_{0,0}$
+# multiplied by `√π` is:
+U[1, 1]*sqrt(π)
+
+# Using an orthonormal basis, the integral of $[f(x,y)]^2$ over the disk is
+# approximately the square of the 2-norm of the coefficients:
+norm(U)^2 - π/(2*sqrt(2))*log1p(sqrt(2))

src/FastTransforms.jl

@@ -35,15 +35,16 @@ import LinearAlgebra: mul!, lmul!, ldiv!
 export leg2cheb, cheb2leg, ultra2ultra, jac2jac,
        lag2lag, jac2ultra, ultra2jac, jac2cheb,
        cheb2jac, ultra2cheb, cheb2ultra,
-       sph2fourier, sphv2fourier, disk2cxf, tri2cheb, tet2cheb,
-       fourier2sph, fourier2sphv, cxf2disk, cheb2tri, cheb2tet
+       sph2fourier, sphv2fourier, disk2cxf, rectdisk2cheb, tri2cheb, tet2cheb,
+       fourier2sph, fourier2sphv, cxf2disk, cheb2rectdisk, cheb2tri, cheb2tet
 
 export plan_leg2cheb, plan_cheb2leg, plan_ultra2ultra, plan_jac2jac,
        plan_lag2lag, plan_jac2ultra, plan_ultra2jac, plan_jac2cheb,
        plan_cheb2jac, plan_ultra2cheb, plan_cheb2ultra,
        plan_sph2fourier, plan_sph_synthesis, plan_sph_analysis,
        plan_sphv2fourier, plan_sphv_synthesis, plan_sphv_analysis,
        plan_disk2cxf, plan_disk_synthesis, plan_disk_analysis,
+       plan_rectdisk2cheb, plan_rectdisk_synthesis, plan_rectdisk_analysis,
        plan_tri2cheb, plan_tri_synthesis, plan_tri_analysis,
        plan_tet2cheb, plan_tet_synthesis, plan_tet_analysis,
        plan_spinsph2fourier, plan_spinsph_synthesis, plan_spinsph_analysis
@@ -91,6 +92,7 @@ include("gaunt.jl")
 export sphones, sphzeros, sphrand, sphrandn, sphevaluate,
        sphvones, sphvzeros, sphvrand, sphvrandn,
        diskones, diskzeros, diskrand, diskrandn,
+       rectdiskones, rectdiskzeros, rectdiskrand, rectdiskrandn,
        triones, trizeros, trirand, trirandn, trievaluate,
        tetones, tetzeros, tetrand, tetrandn,
        spinsphones, spinsphzeros, spinsphrand, spinsphrandn

src/libfasttransforms.jl

@@ -119,22 +119,25 @@ const CHEB2ULTRA           = 10
 const SPHERE               = 11
 const SPHEREV              = 12
 const DISK                 = 13
-const TRIANGLE             = 14
-const TETRAHEDRON          = 15
-const SPINSPHERE           = 16
-const SPHERESYNTHESIS      = 17
-const SPHEREANALYSIS       = 18
-const SPHEREVSYNTHESIS     = 19
-const SPHEREVANALYSIS      = 20
-const DISKSYNTHESIS        = 21
-const DISKANALYSIS         = 22
-const TRIANGLESYNTHESIS    = 23
-const TRIANGLEANALYSIS     = 24
-const TETRAHEDRONSYNTHESIS = 25
-const TETRAHEDRONANALYSIS  = 26
-const SPINSPHERESYNTHESIS  = 27
-const SPINSPHEREANALYSIS   = 28
-const SPHERICALISOMETRY    = 29
+const RECTDISK             = 14
+const TRIANGLE             = 15
+const TETRAHEDRON          = 16
+const SPINSPHERE           = 17
+const SPHERESYNTHESIS      = 18
+const SPHEREANALYSIS       = 19
+const SPHEREVSYNTHESIS     = 20
+const SPHEREVANALYSIS      = 21
+const DISKSYNTHESIS        = 22
+const DISKANALYSIS         = 23
+const RECTDISKSYNTHESIS    = 24
+const RECTDISKANALYSIS     = 25
+const TRIANGLESYNTHESIS    = 26
+const TRIANGLEANALYSIS     = 27
+const TETRAHEDRONSYNTHESIS = 28
+const TETRAHEDRONANALYSIS  = 29
+const SPINSPHERESYNTHESIS  = 30
+const SPINSPHEREANALYSIS   = 31
+const SPHERICALISOMETRY    = 32
 
 
 let k2s = Dict(LEG2CHEB             => "Legendre--Chebyshev",
@@ -151,6 +154,7 @@ let k2s = Dict(LEG2CHEB             => "Legendre--Chebyshev",
                SPHERE               => "Spherical harmonic--Fourier",
                SPHEREV              => "Spherical vector field--Fourier",
                DISK                 => "Zernike--Chebyshev×Fourier",
+               RECTDISK             => "Dunkl-Xu--Chebyshev²",
                TRIANGLE             => "Proriol--Chebyshev²",
                TETRAHEDRON          => "Proriol--Chebyshev³",
                SPINSPHERE           => "Spin-weighted spherical harmonic--Fourier",
@@ -160,6 +164,8 @@ let k2s = Dict(LEG2CHEB             => "Legendre--Chebyshev",
                SPHEREVANALYSIS      => "FFTW Fourier analysis on the sphere (vector field)",
                DISKSYNTHESIS        => "FFTW Chebyshev×Fourier synthesis on the disk",
                DISKANALYSIS         => "FFTW Chebyshev×Fourier analysis on the disk",
+               RECTDISKSYNTHESIS    => "FFTW Chebyshev synthesis on the rectangularized disk",
+               RECTDISKANALYSIS     => "FFTW Chebyshev analysis on the rectangularized disk",
                TRIANGLESYNTHESIS    => "FFTW Chebyshev synthesis on the triangle",
                TRIANGLEANALYSIS     => "FFTW Chebyshev analysis on the triangle",
                TETRAHEDRONSYNTHESIS => "FFTW Chebyshev synthesis on the tetrahedron",
@@ -201,6 +207,7 @@ show(io::IO, p::FTPlan{T, 1, K}) where {T, K} = print(io, "FastTransforms ", kin
 show(io::IO, p::FTPlan{T, 2, SPHERE}) where T = print(io, "FastTransforms ", kind2string(SPHERE), " plan for $(p.n)×$(2p.n-1)-element array of ", T)
 show(io::IO, p::FTPlan{T, 2, SPHEREV}) where T = print(io, "FastTransforms ", kind2string(SPHEREV), " plan for $(p.n)×$(2p.n-1)-element array of ", T)
 show(io::IO, p::FTPlan{T, 2, DISK}) where T = print(io, "FastTransforms ", kind2string(DISK), " plan for $(p.n)×$(4p.n-3)-element array of ", T)
+show(io::IO, p::FTPlan{T, 2, RECTDISK}) where T = print(io, "FastTransforms ", kind2string(RECTDISK), " plan for $(p.n)×$(p.n)-element array of ", T)
 show(io::IO, p::FTPlan{T, 2, TRIANGLE}) where T = print(io, "FastTransforms ", kind2string(TRIANGLE), " plan for $(p.n)×$(p.n)-element array of ", T)
 show(io::IO, p::FTPlan{T, 3, TETRAHEDRON}) where T = print(io, "FastTransforms ", kind2string(TETRAHEDRON), " plan for $(p.n)×$(p.n)×$(p.n)-element array of ", T)
 show(io::IO, p::FTPlan{T, 2, SPINSPHERE}) where T = print(io, "FastTransforms ", kind2string(SPINSPHERE), " plan for $(p.n)×$(2p.n-1)-element array of ", T)
@@ -246,6 +253,8 @@ destroy_plan(p::FTPlan{Float64, 2, SPHEREVSYNTHESIS}) = ccall((:ft_destroy_spher
 destroy_plan(p::FTPlan{Float64, 2, SPHEREVANALYSIS}) = ccall((:ft_destroy_sphere_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
 destroy_plan(p::FTPlan{Float64, 2, DISKSYNTHESIS}) = ccall((:ft_destroy_disk_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
 destroy_plan(p::FTPlan{Float64, 2, DISKANALYSIS}) = ccall((:ft_destroy_disk_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
+destroy_plan(p::FTPlan{Float64, 2, RECTDISKSYNTHESIS}) = ccall((:ft_destroy_rectdisk_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
+destroy_plan(p::FTPlan{Float64, 2, RECTDISKANALYSIS}) = ccall((:ft_destroy_rectdisk_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
 destroy_plan(p::FTPlan{Float64, 2, TRIANGLESYNTHESIS}) = ccall((:ft_destroy_triangle_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
 destroy_plan(p::FTPlan{Float64, 2, TRIANGLEANALYSIS}) = ccall((:ft_destroy_triangle_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
 destroy_plan(p::FTPlan{Float64, 3, TETRAHEDRONSYNTHESIS}) = ccall((:ft_destroy_tetrahedron_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
@@ -299,7 +308,8 @@ unsafe_convert(::Type{Ptr{mpfr_t}}, p::TransposeFTPlan{T, FTPlan{T, N, K}}) wher
 for f in (:leg2cheb, :cheb2leg, :ultra2ultra, :jac2jac,
           :lag2lag, :jac2ultra, :ultra2jac, :jac2cheb,
           :cheb2jac, :ultra2cheb, :cheb2ultra,
-          :sph2fourier, :sphv2fourier, :disk2cxf, :tri2cheb, :tet2cheb)
+          :sph2fourier, :sphv2fourier, :disk2cxf,
+          :rectdisk2cheb, :tri2cheb, :tet2cheb)
     plan_f = Symbol("plan_", f)
     @eval begin
         $plan_f(x::AbstractArray{T}, y...; z...) where T = $plan_f(T, size(x, 1), y...; z...)
@@ -309,8 +319,8 @@ for f in (:leg2cheb, :cheb2leg, :ultra2ultra, :jac2jac,
 end
 
 for (f, plan_f) in ((:fourier2sph, :plan_sph2fourier), (:fourier2sphv, :plan_sphv2fourier),
-                    (:cxf2disk2, :plan_disk2cxf), (:cheb2tri, :plan_tri2cheb),
-                    (:cheb2tet, :plan_tet2cheb))
+                    (:cxf2disk, :plan_disk2cxf), (:cheb2rectdisk, :plan_rectdisk2cheb),
+                    (:cheb2tri, :plan_tri2cheb), (:cheb2tet, :plan_tet2cheb))
     @eval begin
         $f(x::AbstractArray, y...; z...) = $plan_f(x, y...; z...)\x
     end
@@ -498,11 +508,16 @@ function plan_sphv2fourier(::Type{Float64}, n::Integer)
     return FTPlan{Float64, 2, SPHEREV}(plan, n)
 end
 
-function plan_disk2cxf(::Type{Float64}, n::Integer)
-    plan = ccall((:ft_plan_disk2cxf, libfasttransforms), Ptr{ft_plan_struct}, (Cint, ), n)
+function plan_disk2cxf(::Type{Float64}, n::Integer, α, β)
+    plan = ccall((:ft_plan_disk2cxf, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Float64, Float64), n, α, β)
     return FTPlan{Float64, 2, DISK}(plan, n)
 end
 
+function plan_rectdisk2cheb(::Type{Float64}, n::Integer, β)
+    plan = ccall((:ft_plan_rectdisk2cheb, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Float64), n, β)
+    return FTPlan{Float64, 2, RECTDISK}(plan, n)
+end
+
 function plan_tri2cheb(::Type{Float64}, n::Integer, α, β, γ)
     plan = ccall((:ft_plan_tri2cheb, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Float64, Float64, Float64), n, α, β, γ)
     return FTPlan{Float64, 2, TRIANGLE}(plan, n)
@@ -524,6 +539,8 @@ for (fJ, fC, fE, K) in ((:plan_sph_synthesis, :ft_plan_sph_synthesis, :ft_execut
                     (:plan_sphv_analysis, :ft_plan_sphv_analysis, :ft_execute_sphv_analysis, SPHEREVANALYSIS),
                     (:plan_disk_synthesis, :ft_plan_disk_synthesis, :ft_execute_disk_synthesis, DISKSYNTHESIS),
                     (:plan_disk_analysis, :ft_plan_disk_analysis, :ft_execute_disk_analysis, DISKANALYSIS),
+                    (:plan_rectdisk_synthesis, :ft_plan_rectdisk_synthesis, :ft_execute_rectdisk_synthesis, RECTDISKSYNTHESIS),
+                    (:plan_rectdisk_analysis, :ft_plan_rectdisk_analysis, :ft_execute_rectdisk_analysis, RECTDISKANALYSIS),
                     (:plan_tri_synthesis, :ft_plan_tri_synthesis, :ft_execute_tri_synthesis, TRIANGLESYNTHESIS),
                     (:plan_tri_analysis, :ft_plan_tri_analysis, :ft_execute_tri_analysis, TRIANGLEANALYSIS))
     @eval begin
@@ -718,6 +735,8 @@ for (fJ, fC, K) in ((:lmul!, :ft_execute_sph2fourier, SPHERE),
                     (:ldiv!, :ft_execute_fourier2sphv, SPHEREV),
                     (:lmul!, :ft_execute_disk2cxf, DISK),
                     (:ldiv!, :ft_execute_cxf2disk, DISK),
+                    (:lmul!, :ft_execute_rectdisk2cheb, RECTDISK),
+                    (:ldiv!, :ft_execute_cheb2rectdisk, RECTDISK),
                     (:lmul!, :ft_execute_tri2cheb, TRIANGLE),
                     (:ldiv!, :ft_execute_cheb2tri, TRIANGLE))
     @eval begin

src/specialfunctions.jl

@@ -551,6 +551,11 @@ end
 
 trizeros(::Type{T}, m::Int, n::Int) where T = zeros(T, m, n)
 
+const rectdiskrand = trirand
+const rectdiskrandn = trirandn
+const rectdiskones = triones
+const rectdiskzeros = trizeros
+
 """
 Pointwise evaluation of triangular harmonic:
 

test/libfasttransformstests.jl

@@ -105,16 +105,27 @@ FastTransforms.set_num_threads(ceil(Int, Base.Sys.CPU_THREADS/2))
     test_nd_plans(p, ps, pa, A)
 
     A = diskones(Float64, n, 4n-3)
-    p = plan_disk2cxf(A)
+    p = plan_disk2cxf(A, α, β)
     ps = plan_disk_synthesis(A)
     pa = plan_disk_analysis(A)
     test_nd_plans(p, ps, pa, A)
     A = diskones(Float64, n, 4n-3) + im*diskones(Float64, n, 4n-3)
-    p = plan_disk2cxf(A)
+    p = plan_disk2cxf(A, α, β)
     ps = plan_disk_synthesis(A)
     pa = plan_disk_analysis(A)
     test_nd_plans(p, ps, pa, A)
 
+    A = rectdiskones(Float64, n, n)
+    p = plan_rectdisk2cheb(A, β)
+    ps = plan_rectdisk_synthesis(A)
+    pa = plan_rectdisk_analysis(A)
+    test_nd_plans(p, ps, pa, A)
+    A = rectdiskones(Float64, n, n) + im*rectdiskones(Float64, n, n)
+    p = plan_rectdisk2cheb(A, β)
+    ps = plan_rectdisk_synthesis(A)
+    pa = plan_rectdisk_analysis(A)
+    test_nd_plans(p, ps, pa, A)
+
     A = triones(Float64, n, n)
     p = plan_tri2cheb(A, α, β, γ)
     ps = plan_tri_synthesis(A)