let the build script point to master

add wrapper and examples

Author: MikaelSlevinsky

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/sphericalisometries.jl

@@ -0,0 +1,115 @@
+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_ZY_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)
+
+
+# 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
 

src/libfasttransforms.jl

@@ -134,6 +134,7 @@ const TETRAHEDRONSYNTHESIS = 25
 const TETRAHEDRONANALYSIS  = 26
 const SPINSPHERESYNTHESIS  = 27
 const SPINSPHEREANALYSIS   = 28
+const SPHERICALISOMETRY    = 29
 
 
 let k2s = Dict(LEG2CHEB             => "Legendre--Chebyshev",
@@ -164,7 +165,8 @@ let k2s = Dict(LEG2CHEB             => "Legendre--Chebyshev",
                TETRAHEDRONSYNTHESIS => "FFTW Chebyshev synthesis 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)")
+               SPINSPHEREANALYSIS   => "FFTW Fourier analysis on the sphere (spin-weighted)",
+               SPHERICALISOMETRY    => "Spherical isometry")
     global kind2string
     kind2string(k::Integer) = k2s[Int(k)]
 end
@@ -204,6 +206,7 @@ show(io::IO, p::FTPlan{T, 3, TETRAHEDRON}) where T = print(io, "FastTransforms "
 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)
+show(io::IO, p::FTPlan{T, 2, SPHERICALISOMETRY}) where T = print(io, "FastTransforms ", kind2string(SPHERICALISOMETRY), " plan for $(p.n)×$(2p.n-1)-element array of ", T)
 
 function checksize(p::FTPlan{T}, x::Array{T}) where T
     if p.n != size(x, 1)
@@ -222,6 +225,12 @@ for K in (SPHERE, SPHEREV, DISK, SPINSPHERE)
     end
 end
 
+function checksize(p::FTPlan{T, 2, SPHERICALISOMETRY}, x::Matrix{T}) where T
+    if p.n != size(x, 1) || 2p.n-1 != size(x, 2)
+        throw(DimensionMismatch("This FTPlan must operate on arrays of size $(p.n) × $(2p.n-1)."))
+    end
+end
+
 unsafe_convert(::Type{Ptr{ft_plan_struct}}, p::FTPlan) = p.plan
 unsafe_convert(::Type{Ptr{mpfr_t}}, p::FTPlan) = unsafe_convert(Ptr{mpfr_t}, p.plan)
 
@@ -243,6 +252,7 @@ destroy_plan(p::FTPlan{Float64, 3, TETRAHEDRONSYNTHESIS}) = ccall((:ft_destroy_t
 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)
+destroy_plan(p::FTPlan{Float64, 2, SPHERICALISOMETRY}) = ccall((:ft_destroy_sph_isometry_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
 
 struct AdjointFTPlan{T, S}
     parent::S
@@ -595,6 +605,11 @@ function lmul!(p::FTPlan{Complex{Float64}, 2, SPINSPHEREANALYSIS}, x::Matrix{Com
     return x
 end
 
+function plan_sph_isometry(::Type{Float64}, n::Integer)
+    plan = ccall((:ft_plan_sph_isometry, libfasttransforms), Ptr{ft_plan_struct}, (Cint, ), n)
+    return FTPlan{Float64, 2, SPHERICALISOMETRY}(plan, n)
+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))
 *(p::TransposeFTPlan{T}, x::Array{T}) where T = lmul!(p, deepcopy(x))
@@ -738,6 +753,59 @@ function ldiv!(p::FTPlan{Complex{Float64}, 2, SPINSPHERE}, x::Matrix{Complex{Flo
     return x
 end
 
+function execute_sph_polar_rotation!(x::Matrix{Float64}, α)
+    ccall((:ft_execute_sph_polar_rotation, libfasttransforms), Cvoid, (Ptr{Float64}, Cint, Cint, Float64, Float64), x, size(x, 1), size(x, 2), sin(α), cos(α))
+    return x
+end
+
+function execute_sph_polar_reflection!(x::Matrix{Float64})
+    ccall((:ft_execute_sph_polar_reflection, libfasttransforms), Cvoid, (Ptr{Float64}, Cint, Cint), x, size(x, 1), size(x, 2))
+    return x
+end
+
+struct ft_orthogonal_transformation
+    Q::NTuple{9, Float64}
+end
+
+function convert(::Type{ft_orthogonal_transformation}, Q::AbstractMatrix)
+    @assert size(Q, 1) ≥ 3 && size(Q, 2) ≥ 3
+    return ft_orthogonal_transformation((Q[1, 1], Q[2, 1], Q[3, 1], Q[1, 2], Q[2, 2], Q[3, 2], Q[1, 3], Q[2, 3], Q[3, 3]))
+end
+
+function execute_sph_orthogonal_transformation!(p::FTPlan{Float64, 2, SPHERICALISOMETRY}, Q, x::Matrix{Float64})
+    checksize(p, x)
+    ccall((:ft_execute_sph_orthogonal_transformation, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ft_orthogonal_transformation, Ptr{Float64}, Cint, Cint), p, Q, x, size(x, 1), size(x, 2))
+    return x
+end
+
+function execute_sph_ZY_axis_exchange!(p::FTPlan{Float64, 2, SPHERICALISOMETRY}, x::Matrix{Float64})
+    checksize(p, x)
+    ccall((:ft_execute_sph_ZY_axis_exchange, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, Ptr{Float64}, Cint, Cint), p, x, size(x, 1), size(x, 2))
+    return x
+end
+
+function execute_sph_rotation!(p::FTPlan{Float64, 2, SPHERICALISOMETRY}, α, β, γ, x::Matrix{Float64})
+    checksize(p, x)
+    ccall((:ft_execute_sph_rotation, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, Float64, Float64, Float64, Ptr{Float64}, Cint, Cint), p, α, β, γ, x, size(x, 1), size(x, 2))
+    return x
+end
+
+struct ft_reflection
+    w::NTuple{3, Float64}
+end
+
+function convert(::Type{ft_reflection}, w::AbstractVector)
+    @assert length(w) ≥ 3
+    return ft_reflection((w[1], w[2], w[3]))
+end
+
+function execute_sph_reflection!(p::FTPlan{Float64, 2, SPHERICALISOMETRY}, w, x::Matrix{Float64})
+    checksize(p, x)
+    ccall((:ft_execute_sph_reflection, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ft_reflection, Ptr{Float64}, Cint, Cint), p, w, x, size(x, 1), size(x, 2))
+    return x
+end
+execute_sph_reflection!(p::FTPlan{Float64, 2, SPHERICALISOMETRY}, w1, w2, w3, x::Matrix{Float64}) = execute_sph_reflection!(p, ft_reflection(w1, w2, w3), x)
+
 *(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))