[FastTransforms] v0.4 (#125)

* let the build script point to master

add wrapper and examples

* add extra conversions from NTuple{*, Float64}

* rename ZY_ to yz_

* update disk2cxf, add rectdisk2cheb

* remove nightly tests

* remove unnecessary parentheses

Author: MikaelSlevinsky

.github/workflows/ci.yml

@@ -11,7 +11,6 @@ jobs:
         version:
           - '1.3'
           - '1.5'
-          - 'nightly'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -1,6 +1,6 @@
 name = "FastTransforms"
 uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c"
-version = "0.10.3"
+version = "0.11.0"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
@@ -19,15 +19,15 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 
 [compat]
-AbstractFFTs = "0.4, 0.5"
-ArrayLayouts = "0.3.7, 0.4"
+AbstractFFTs = "0.5"
+ArrayLayouts = "0.4"
 BinaryProvider = "0.5"
 DSP = "0.6"
 FFTW = "1"
 FastGaussQuadrature = "0.4"
-FastTransforms_jll = "0.3.3"
-FillArrays = "0.8, 0.9, 0.10"
+FastTransforms_jll = "0.4.0"
+FillArrays = "0.9, 0.10"
 Reexport = "0.2"
-SpecialFunctions = "0.8, 0.9, 0.10"
+SpecialFunctions = "0.10, 1"
 ToeplitzMatrices = "0.6"
 julia = "1.3"

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
 

deps/build.jl

@@ -1,8 +1,6 @@
 using BinaryProvider
 import Libdl
 
-version = v"0.3.3"
-
 const extension = Sys.isapple() ? "dylib" : Sys.islinux() ? "so" : Sys.iswindows() ? "dll" : ""
 
 print_error() = error(
@@ -26,10 +24,11 @@ if ft_build_from_source == "true"
         if [ -d "FastTransforms" ]; then
             cd FastTransforms
             git fetch
-            git checkout v$version
+            git checkout master
+            git pull
             cd ..
         else
-            git clone -b v$version https://github.com/MikaelSlevinsky/FastTransforms.git FastTransforms
+            git clone https://github.com/MikaelSlevinsky/FastTransforms.git FastTransforms
         fi
         cd FastTransforms
         $make assembly

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))

examples/sphericalisometries.jl

@@ -0,0 +1,124 @@
+function threshold!(A::AbstractArray, ϵ)
+    for i in eachindex(A)
+        if abs(A[i]) < ϵ A[i] = 0 end
+    end
+    A
+end
+
+using FastTransforms, LinearAlgebra, Random, Test
+
+# The colatitudinal grid (mod π):
+N = 10
+θ = (0.5:N-0.5)/N
+
+# The longitudinal grid (mod π):
+M = 2*N-1
+φ = (0:M-1)*2/M
+
+x = [cospi(φ)*sinpi(θ) for θ in θ, φ in φ]
+y = [sinpi(φ)*sinpi(θ) for θ in θ, φ in φ]
+z = [cospi(θ) for θ in θ, φ in φ]
+
+P = plan_sph2fourier(Float64, N)
+PA = plan_sph_analysis(Float64, N, M)
+J = FastTransforms.plan_sph_isometry(Float64, N)
+
+
+f = (x, y, z) -> x^2+y^4+x^2*y*z^3-x*y*z^2
+
+
+F = f.(x, y, z)
+V = PA*F
+U = threshold!(P\V, 100eps())
+FastTransforms.execute_sph_yz_axis_exchange!(J, U)
+FR = f.(x, -z, -y)
+VR = PA*FR
+UR = threshold!(P\VR, 100eps())
+@test U ≈ UR
+norm(U-UR)
+
+
+α, β, γ = 0.123, 0.456, 0.789
+
+# Isometry built up from ZYZR
+A = [cos(α) -sin(α) 0; sin(α) cos(α) 0; 0 0 1]
+B = [cos(β) 0 -sin(β); 0 1 0; sin(β) 0 cos(β)]
+C = [cos(γ) -sin(γ) 0; sin(γ) cos(γ) 0; 0 0 1]
+R = diagm([1, 1, 1.0])
+Q = A*B*C*R
+
+# Transform the sampling grid. Note that `Q` is transposed here.
+u = Q[1,1]*x + Q[2,1]*y + Q[3,1]*z
+v = Q[1,2]*x + Q[2,2]*y + Q[3,2]*z
+w = Q[1,3]*x + Q[2,3]*y + Q[3,3]*z
+
+F = f.(x, y, z)
+V = PA*F
+U = threshold!(P\V, 100eps())
+FastTransforms.execute_sph_rotation!(J, α, β, γ, U)
+FR = f.(u, v, w)
+VR = PA*FR
+UR = threshold!(P\VR, 100eps())
+@test U ≈ UR
+norm(U-UR)
+
+
+F = f.(x, y, z)
+V = PA*F
+U = threshold!(P\V, 100eps())
+FastTransforms.execute_sph_polar_reflection!(U)
+FR = f.(x, y, -z)
+VR = PA*FR
+UR = threshold!(P\VR, 100eps())
+@test U ≈ UR
+norm(U-UR)
+
+
+# Isometry built up from planar reflection
+W = [0.123, 0.456, 0.789]
+H = w -> I - 2/(w'w)*w*w'
+Q = H(W)
+
+# Transform the sampling grid. Note that `Q` is transposed here.
+u = Q[1,1]*x + Q[2,1]*y + Q[3,1]*z
+v = Q[1,2]*x + Q[2,2]*y + Q[3,2]*z
+w = Q[1,3]*x + Q[2,3]*y + Q[3,3]*z
+
+F = f.(x, y, z)
+V = PA*F
+U = threshold!(P\V, 100eps())
+FastTransforms.execute_sph_reflection!(J, W, U)
+FR = f.(u, v, w)
+VR = PA*FR
+UR = threshold!(P\VR, 100eps())
+@test U ≈ UR
+norm(U-UR)
+
+F = f.(x, y, z)
+V = PA*F
+U = threshold!(P\V, 100eps())
+FastTransforms.execute_sph_reflection!(J, (W[1], W[2], W[3]), U)
+FR = f.(u, v, w)
+VR = PA*FR
+UR = threshold!(P\VR, 100eps())
+@test U ≈ UR
+norm(U-UR)
+
+# Random orthogonal transformation
+Random.seed!(0)
+Q = qr(rand(3, 3)).Q
+
+# Transform the sampling grid, note that `Q` is transposed here.
+u = Q[1,1]*x + Q[2,1]*y + Q[3,1]*z
+v = Q[1,2]*x + Q[2,2]*y + Q[3,2]*z
+w = Q[1,3]*x + Q[2,3]*y + Q[3,3]*z
+
+F = f.(x, y, z)
+V = PA*F
+U = threshold!(P\V, 100eps())
+FastTransforms.execute_sph_orthogonal_transformation!(J, Q, U)
+FR = f.(u, v, w)
+VR = PA*FR
+UR = threshold!(P\VR, 100eps())
+@test U ≈ UR
+norm(U-UR)

src/FastTransforms.jl

@@ -9,7 +9,7 @@ import DSP
 @reexport using FFTW
 
 import Base: unsafe_convert, eltype, ndims, adjoint, transpose, show, *, \,
-             inv, length, size, view, getindex
+             inv, length, size, view, getindex, convert
 
 import Base.GMP: Limb
 
@@ -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