Add spin-weighted spherical harmonics (#103)

* add spin-weighted spherical harmonics

* update freebsd image

* add horner and clenshaw routines

don't export

Author: MikaelSlevinsky

.cirrus.yml

@@ -1,5 +1,5 @@
 freebsd_instance:
-  image: freebsd-12-0-release-amd64
+  image: freebsd-12-1-release-amd64
 task:
   name: FreeBSD
   env:

Project.toml

@@ -1,6 +1,6 @@
 name = "FastTransforms"
 uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c"
-version = "0.8.2"
+version = "0.9.0"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
@@ -22,7 +22,7 @@ BinaryProvider = "0.5.8"
 DSP = "0.6"
 FFTW = "1"
 FastGaussQuadrature = "0.4"
-FastTransforms_jll = "0.2.13"
+FastTransforms_jll = "0.3.0"
 Reexport = "0.2"
 SpecialFunctions = "0.8, 0.9, 0.10"
 ToeplitzMatrices = "0.6"

README.md

@@ -19,7 +19,7 @@ julia> using FastTransforms, LinearAlgebra
 
 ## Fast orthogonal polynomial transforms
 
-The 26 orthogonal polynomial transforms are listed in `FastTransforms.kind2string.(0:25)`. 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 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 Chebyshev--Legendre transform
 

examples/spinweighted.jl

@@ -0,0 +1,64 @@
+#############
+# This example plays with analysis of:
+#
+#   f(r) = e^{i k⋅r},
+#
+# for some k ∈ ℝ³ and where r ∈ 𝕊², using spin-0 spherical harmonics.
+#
+# It applies ð, the spin-raising operator,
+# both on the spin-0 coefficients as well as the original function,
+# followed by a spin-1 analysis to compare coefficients.
+#
+# See also sphere.jl
+# For the storage pattern of the arrays, please consult the documentation.
+#############
+
+using FastTransforms, LinearAlgebra
+
+# 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
+
+k = [2/7, 3/7, 6/7]
+r = (θ,φ) -> [sinpi(θ)*cospi(φ), sinpi(θ)*sinpi(φ), cospi(θ)]
+
+# On the tensor product grid, our function samples are:
+
+F = [exp(im*(k⋅r(θ,φ))) for θ in θ, φ in φ]
+
+P = plan_spinsph2fourier(F, 0)
+PA = plan_spinsph_analysis(F, 0)
+
+# Its spin-0 spherical harmonic coefficients are:
+
+U⁰ = P\(PA*F)
+
+norm(U⁰) ≈ sqrt(4π)
+
+# Spin can be incremented by applying ð, either on the spin-0 coefficients:
+
+U¹c = zero(U⁰)
+for n in 1:N-1
+    U¹c[n, 1] = sqrt(n*(n+1))*U⁰[n+1, 1]
+end
+for m in 1:M÷2
+    for n in 0:N-1
+        U¹c[n+1, 2m] = -sqrt((n+m)*(n+m+1))*U⁰[n+1, 2m]
+        U¹c[n+1, 2m+1] = sqrt((n+m)*(n+m+1))*U⁰[n+1, 2m+1]
+    end
+end
+
+# or on the original function through analysis with spin-1 spherical harmonics:
+
+F = [-(k[1]*(im*cospi(θ)*cospi(φ) + sinpi(φ)) + k[2]*(im*cospi(θ)*sinpi(φ)-cospi(φ)) - im*k[3]*sinpi(θ))*exp(im*(k⋅r(θ,φ))) for θ in θ, φ in φ]
+
+P = plan_spinsph2fourier(F, 1)
+PA = plan_spinsph_analysis(F, 1)
+
+U¹s = P\(PA*F)
+
+norm(U¹c) ≈ norm(U¹s) ≈ sqrt(8π/3*(k⋅k))

src/FastTransforms.jl

@@ -44,7 +44,8 @@ export plan_leg2cheb, plan_cheb2leg, plan_ultra2ultra, plan_jac2jac,
        plan_sphv2fourier, plan_sphv_synthesis, plan_sphv_analysis,
        plan_disk2cxf, plan_disk_synthesis, plan_disk_analysis,
        plan_tri2cheb, plan_tri_synthesis, plan_tri_analysis,
-       plan_tet2cheb, plan_tet_synthesis, plan_tet_analysis
+       plan_tet2cheb, plan_tet_synthesis, plan_tet_analysis,
+       plan_spinsph2fourier, plan_spinsph_synthesis, plan_spinsph_analysis
 
 include("libfasttransforms.jl")
 
@@ -87,7 +88,8 @@ export sphones, sphzeros, sphrand, sphrandn, sphevaluate,
        sphvones, sphvzeros, sphvrand, sphvrandn,
        diskones, diskzeros, diskrand, diskrandn,
        triones, trizeros, trirand, trirandn, trievaluate,
-       tetones, tetzeros, tetrand, tetrandn
+       tetones, tetzeros, tetrand, tetrandn,
+       spinsphones, spinsphzeros, spinsphrand, spinsphrandn
 
 lgamma(x) = logabsgamma(x)[1]
 

src/libfasttransforms.jl

@@ -43,6 +43,38 @@ end
 
 set_num_threads(n::Integer) = ccall((:ft_set_num_threads, libfasttransforms), Cvoid, (Cint, ), n)
 
+function horner!(c::Vector{Float64}, x::Vector{Float64}, f::Vector{Float64})
+    @assert length(x) == length(f)
+    ccall((:ft_horner, libfasttransforms), Cvoid, (Cint, Ptr{Float64}, Cint, Cint, Ptr{Float64}, Ptr{Float64}), length(c), c, 1, length(x), x, f)
+end
+
+function horner!(c::Vector{Float32}, x::Vector{Float32}, f::Vector{Float32})
+    @assert length(x) == length(f)
+    ccall((:ft_hornerf, libfasttransforms), Cvoid, (Cint, Ptr{Float32}, Cint, Cint, Ptr{Float32}, Ptr{Float32}), length(c), c, 1, length(x), x, f)
+end
+
+function clenshaw!(c::Vector{Float64}, x::Vector{Float64}, f::Vector{Float64})
+    @assert length(x) == length(f)
+    ccall((:ft_clenshaw, libfasttransforms), Cvoid, (Cint, Ptr{Float64}, Cint, Cint, Ptr{Float64}, Ptr{Float64}), length(c), c, 1, length(x), x, f)
+end
+
+function clenshaw!(c::Vector{Float32}, x::Vector{Float32}, f::Vector{Float32})
+    @assert length(x) == length(f)
+    ccall((:ft_clenshawf, libfasttransforms), Cvoid, (Cint, Ptr{Float32}, Cint, Cint, Ptr{Float32}, Ptr{Float32}), length(c), c, 1, length(x), x, f)
+end
+
+function clenshaw!(c::Vector{Float64}, A::Vector{Float64}, B::Vector{Float64}, C::Vector{Float64}, x::Vector{Float64}, phi0::Vector{Float64}, f::Vector{Float64})
+    @assert length(c) == length(A) == length(B) == length(C)-1
+    @assert length(x) == length(phi0) == length(f)
+    ccall((:ft_orthogonal_polynomial_clenshaw, libfasttransforms), Cvoid, (Cint, Ptr{Float64}, Cint, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Cint, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}), length(c), c, 1, A, B, C, length(x), x, phi0, f)
+end
+
+function clenshaw!(c::Vector{Float32}, A::Vector{Float32}, B::Vector{Float32}, C::Vector{Float32}, x::Vector{Float32}, phi0::Vector{Float32}, f::Vector{Float32})
+    @assert length(c) == length(A) == length(B) == length(C)-1
+    @assert length(x) == length(phi0) == length(f)
+    ccall((:ft_orthogonal_polynomial_clenshawf, libfasttransforms), Cvoid, (Cint, Ptr{Float32}, Cint, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Cint, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}), length(c), c, 1, A, B, C, length(x), x, phi0, f)
+end
+
 const LEG2CHEB              = 0
 const CHEB2LEG              = 1
 const ULTRA2ULTRA           = 2
@@ -59,16 +91,20 @@ const SPHEREV              = 12
 const DISK                 = 13
 const TRIANGLE             = 14
 const TETRAHEDRON          = 15
-const SPHERESYNTHESIS      = 16
-const SPHEREANALYSIS       = 17
-const SPHEREVSYNTHESIS     = 18
-const SPHEREVANALYSIS      = 19
-const DISKSYNTHESIS        = 20
-const DISKANALYSIS         = 21
-const TRIANGLESYNTHESIS    = 22
-const TRIANGLEANALYSIS     = 23
-const TETRAHEDRONSYNTHESIS = 24
-const TETRAHEDRONANALYSIS  = 25
+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
+
 
 let k2s = Dict(LEG2CHEB             => "Legendre--Chebyshev",
                CHEB2LEG             => "Chebyshev--Legendre",
@@ -86,6 +122,7 @@ let k2s = Dict(LEG2CHEB             => "Legendre--Chebyshev",
                DISK                 => "Zernike--Chebyshev×Fourier",
                TRIANGLE             => "Proriol--Chebyshev²",
                TETRAHEDRON          => "Proriol--Chebyshev³",
+               SPINSPHERE           => "Spin-weighted spherical harmonic--Fourier",
                SPHERESYNTHESIS      => "FFTW Fourier synthesis on the sphere",
                SPHEREANALYSIS       => "FFTW Fourier analysis on the sphere",
                SPHEREVSYNTHESIS     => "FFTW Fourier synthesis on the sphere (vector field)",
@@ -95,7 +132,9 @@ let k2s = Dict(LEG2CHEB             => "Legendre--Chebyshev",
                TRIANGLESYNTHESIS    => "FFTW Chebyshev synthesis on the triangle",
                TRIANGLEANALYSIS     => "FFTW Chebyshev analysis on the triangle",
                TETRAHEDRONSYNTHESIS => "FFTW Chebyshev synthesis on the tetrahedron",
-               TETRAHEDRONANALYSIS  => "FFTW Chebyshev analysis on the tetrahedron")
+               TETRAHEDRONANALYSIS  => "FFTW Chebyshev analysis on the tetrahedron",
+               SPINSPHERESYNTHESIS  => "FFTW Fourier synthesis on the sphere (spin-weighted)",
+               SPINSPHEREANALYSIS   => "FFTW Fourier analysis on the sphere (spin-weighted)")
     global kind2string
     kind2string(k::Integer) = k2s[Int(k)]
 end
@@ -132,6 +171,7 @@ show(io::IO, p::FTPlan{T, 2, SPHEREV}) where T = print(io, "FastTransforms ", ki
 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, 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)
 show(io::IO, p::FTPlan{T, 2, K}) where {T, K} = print(io, "FastTransforms plan for ", kind2string(K), " for $(p.n)×$(p.m)-element array of ", T)
 show(io::IO, p::FTPlan{T, 3, K}) where {T, K} = print(io, "FastTransforms plan for ", kind2string(K), " for $(p.n)×$(p.l)×$(p.m)-element array of ", T)
 
@@ -141,7 +181,7 @@ function checksize(p::FTPlan{T}, x::Array{T}) where T
     end
 end
 
-for K in (SPHERE, SPHEREV, DISK)
+for K in (SPHERE, SPHEREV, DISK, SPINSPHERE)
     @eval function checksize(p::FTPlan{T, 2, $K}, x::Matrix{T}) where T
         if p.n != size(x, 1)
             throw(DimensionMismatch("FTPlan has dimensions $(p.n) × $(p.n), x has leading dimension $(size(x, 1))"))
@@ -160,6 +200,7 @@ destroy_plan(p::FTPlan{Float64, 1}) = ccall((:ft_destroy_tb_eigen_FMM, libfasttr
 destroy_plan(p::FTPlan{BigFloat, 1}) = ccall((:ft_mpfr_destroy_plan, libfasttransforms), Cvoid, (Ptr{mpfr_t}, Cint), p, p.n)
 destroy_plan(p::FTPlan{Float64, 2}) = ccall((:ft_destroy_harmonic_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
 destroy_plan(p::FTPlan{Float64, 3}) = ccall((:ft_destroy_tetrahedral_harmonic_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
+destroy_plan(p::FTPlan{Complex{Float64}, 2, SPINSPHERE}) = ccall((:ft_destroy_spin_harmonic_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
 destroy_plan(p::FTPlan{Float64, 2, SPHERESYNTHESIS}) = ccall((:ft_destroy_sphere_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
 destroy_plan(p::FTPlan{Float64, 2, SPHEREANALYSIS}) = ccall((:ft_destroy_sphere_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
 destroy_plan(p::FTPlan{Float64, 2, SPHEREVSYNTHESIS}) = ccall((:ft_destroy_sphere_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
@@ -170,6 +211,8 @@ destroy_plan(p::FTPlan{Float64, 2, TRIANGLESYNTHESIS}) = ccall((:ft_destroy_tria
 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)
 destroy_plan(p::FTPlan{Float64, 3, TETRAHEDRONANALYSIS}) = ccall((:ft_destroy_tetrahedron_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
+destroy_plan(p::FTPlan{Complex{Float64}, 2, SPINSPHERESYNTHESIS}) = ccall((:ft_destroy_spinsphere_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
+destroy_plan(p::FTPlan{Complex{Float64}, 2, SPINSPHEREANALYSIS}) = ccall((:ft_destroy_spinsphere_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
 
 struct AdjointFTPlan{T, S}
     parent::S
@@ -233,6 +276,9 @@ for (f, plan_f) in ((:fourier2sph, :plan_sph2fourier), (:fourier2sphv, :plan_sph
     end
 end
 
+plan_spinsph2fourier(x::AbstractArray{T}, y...; z...) where T = plan_spinsph2fourier(T, size(x, 1), y...; z...)
+spinsph2fourier(x::AbstractArray, y...; z...) = plan_spinsph2fourier(x, y...; z...)*x
+fourier2spinsph(x::AbstractArray, y...; z...) = plan_spinsph2fourier(x, y...; z...)\x
 
 function plan_leg2cheb(::Type{Float32}, n::Integer; normleg::Bool=false, normcheb::Bool=false)
     plan = ccall((:ft_plan_legendre_to_chebyshevf, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint), normleg, normcheb, n)
@@ -427,6 +473,11 @@ function plan_tet2cheb(::Type{Float64}, n::Integer, α, β, γ, δ)
     return FTPlan{Float64, 3, TETRAHEDRON}(plan, n)
 end
 
+function plan_spinsph2fourier(::Type{Complex{Float64}}, n::Integer, s::Integer)
+    plan = ccall((:ft_plan_spinsph2fourier, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint), n, s)
+    return FTPlan{Complex{Float64}, 2, SPINSPHERE}(plan, n)
+end
+
 for (fJ, fC, fE, K) in ((:plan_sph_synthesis, :ft_plan_sph_synthesis, :ft_execute_sph_synthesis, SPHERESYNTHESIS),
                     (:plan_sph_analysis, :ft_plan_sph_analysis, :ft_execute_sph_analysis, SPHEREANALYSIS),
                     (:plan_sphv_synthesis, :ft_plan_sphv_synthesis, :ft_execute_sphv_synthesis, SPHEREVSYNTHESIS),
@@ -484,6 +535,35 @@ function lmul!(p::FTPlan{Float64, 3, TETRAHEDRONANALYSIS}, x::Array{Float64, 3})
     return x
 end
 
+plan_spinsph_synthesis(x::Matrix{T}, s::Integer) where T = plan_spinsph_synthesis(T, size(x, 1), size(x, 2), s)
+
+function plan_spinsph_synthesis(::Type{Complex{Float64}}, n::Integer, m::Integer, s::Integer)
+    plan = ccall((:ft_plan_spinsph_synthesis, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint), n, m, s)
+    return FTPlan{Complex{Float64}, 2, SPINSPHERESYNTHESIS}(plan, n, m)
+end
+
+function lmul!(p::FTPlan{Complex{Float64}, 2, SPINSPHERESYNTHESIS}, x::Matrix{Complex{Float64}})
+    if p.n != size(x, 1) || p.m != size(x, 2)
+        throw(DimensionMismatch("FTPlan has dimensions $(p.n) × $(p.m), x has dimensions $(size(x, 1)) × $(size(x, 2))"))
+    end
+    ccall((:ft_execute_spinsph_synthesis, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, Ptr{Float64}, Cint, Cint), p, x, size(x, 1), size(x, 2))
+    return x
+end
+
+plan_spinsph_analysis(x::Matrix{T}, s::Integer) where T = plan_spinsph_analysis(T, size(x, 1), size(x, 2), s)
+
+function plan_spinsph_analysis(::Type{Complex{Float64}}, n::Integer, m::Integer, s::Integer)
+    plan = ccall((:ft_plan_spinsph_analysis, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint), n, m, s)
+    return FTPlan{Complex{Float64}, 2, SPINSPHEREANALYSIS}(plan, n, m)
+end
+
+function lmul!(p::FTPlan{Complex{Float64}, 2, SPINSPHEREANALYSIS}, x::Matrix{Complex{Float64}})
+    if p.n != size(x, 1) || p.m != size(x, 2)
+        throw(DimensionMismatch("FTPlan has dimensions $(p.n) × $(p.m), x has dimensions $(size(x, 1)) × $(size(x, 2))"))
+    end
+    ccall((:ft_execute_spinsph_analysis, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, Ptr{Float64}, Cint, Cint), p, x, size(x, 1), size(x, 2))
+    return x
+end
 
 *(p::FTPlan{T}, x::Array{T}) where T = lmul!(p, deepcopy(x))
 *(p::AdjointFTPlan{T}, x::Array{T}) where T = lmul!(p, deepcopy(x))
@@ -634,6 +714,18 @@ function ldiv!(p::FTPlan{Float64, 3, TETRAHEDRON}, x::Array{Float64, 3})
     return x
 end
 
+function lmul!(p::FTPlan{Complex{Float64}, 2, SPINSPHERE}, x::Matrix{Complex{Float64}})
+    checksize(p, x)
+    ccall((:ft_execute_spinsph2fourier, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, Ptr{Complex{Float64}}, Cint, Cint), p, x, size(x, 1), size(x, 2))
+    return x
+end
+
+function ldiv!(p::FTPlan{Complex{Float64}, 2, SPINSPHERE}, x::Matrix{Complex{Float64}})
+    checksize(p, x)
+    ccall((:ft_execute_fourier2spinsph, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, Ptr{Complex{Float64}}, Cint, Cint), p, x, size(x, 1), size(x, 2))
+    return x
+end
+
 *(p::FTPlan{T}, x::Array{Complex{T}}) where T = lmul!(p, deepcopy(x))
 *(p::AdjointFTPlan{T}, x::Array{Complex{T}}) where T = lmul!(p, deepcopy(x))
 *(p::TransposeFTPlan{T}, x::Array{Complex{T}}) where T = lmul!(p, deepcopy(x))

src/specialfunctions.jl

@@ -605,3 +605,50 @@ function tetones(::Type{T}, l::Int, m::Int, n::Int) where T
 end
 
 tetzeros(::Type{T}, l::Int, m::Int, n::Int) where T = zeros(T, l, m, n)
+
+function spinsphrand(::Type{T}, m::Int, n::Int, s::Int) where T
+    A = zeros(T, m, n)
+    as = abs(s)
+    for i = 1:m-as
+        A[i,1] = rand(T)
+    end
+    for j = 1:n÷2
+        for i = 1:m-max(j, as)
+            A[i,2j] = rand(T)
+            A[i,2j+1] = rand(T)
+        end
+    end
+    A
+end
+
+function spinsphrandn(::Type{T}, m::Int, n::Int, s::Int) where T
+    A = zeros(T, m, n)
+    as = abs(s)
+    for i = 1:m-as
+        A[i,1] = randn(T)
+    end
+    for j = 1:n÷2
+        for i = 1:m-max(j, as)
+            A[i,2j] = randn(T)
+            A[i,2j+1] = randn(T)
+        end
+    end
+    A
+end
+
+function spinsphones(::Type{T}, m::Int, n::Int, s::Int) where T
+    A = zeros(T, m, n)
+    as = abs(s)
+    for i = 1:m-as
+        A[i,1] = one(T)
+    end
+    for j = 1:n÷2
+        for i = 1:m-max(j, as)
+            A[i,2j] = one(T)
+            A[i,2j+1] = one(T)
+        end
+    end
+    A
+end
+
+spinsphzeros(::Type{T}, m::Int, n::Int) where T = zeros(T, m, n)