add first stab at synthesis and analysis

Author: MikaelSlevinsky

src/SphericalHarmonics/SphericalHarmonics.jl

@@ -13,6 +13,7 @@ include("slowplan.jl")
 include("Butterfly.jl")
 include("fastplan.jl")
 include("thinplan.jl")
+include("synthesisanalysis.jl")
 
 function plan_sph2fourier(A::AbstractMatrix; opts...)
     M, N = size(A)

src/SphericalHarmonics/slowplan.jl

@@ -58,7 +58,7 @@ end
 @inline length(P::Pnmp2toPlm) = length(P.rotations)
 @inline getindex(P::Pnmp2toPlm, i::Int) = P.rotations[i]
 
-function Base.A_mul_B!(C::Pnmp2toPlm,A::AbstractVecOrMat)
+function Base.A_mul_B!(C::Pnmp2toPlm, A::AbstractVecOrMat)
     @inbounds for i = 1:length(C)
         A_mul_B!(C.rotations[i], A)
     end

src/SphericalHarmonics/synthesisanalysis.jl

@@ -0,0 +1,190 @@
+immutable SynthesisPlan{T, P1, P2}
+    planθ::P1
+    planφ::P2
+    C::ColumnPermutation
+    temp::Vector{T}
+end
+
+function plan_synthesis{T<:FFTW.fftwNumber}(A::Matrix{T})
+    m, n = size(A)
+    x = FFTW.FakeArray(T, m)
+    y = FFTW.FakeArray(T, n)
+    planθ = FFTW.plan_r2r!(x, FFTW.REDFT01), FFTW.plan_r2r!(x, FFTW.RODFT01)
+    planφ = FFTW.plan_r2r!(y, FFTW.HC2R)
+    C = ColumnPermutation(vcat(1:2:n, 2:2:n))
+    SynthesisPlan(planθ, planφ, C, zeros(T, n))
+end
+
+immutable AnalysisPlan{T, P1, P2}
+    planθ::P1
+    planφ::P2
+    C::ColumnPermutation
+    temp::Vector{T}
+end
+
+function plan_analysis{T<:FFTW.fftwNumber}(A::Matrix{T})
+    m, n = size(A)
+    x = FFTW.FakeArray(T, m)
+    y = FFTW.FakeArray(T, n)
+    planθ = FFTW.plan_r2r!(x, FFTW.REDFT10), FFTW.plan_r2r!(x, FFTW.RODFT10)
+    planφ = FFTW.plan_r2r!(y, FFTW.R2HC)
+    C = ColumnPermutation(vcat(1:2:n, 2:2:n))
+    AnalysisPlan(planθ, planφ, C, zeros(T, n))
+end
+
+function Base.A_mul_B!{T}(Y::Matrix{T}, P::SynthesisPlan{T}, X::Matrix{T})
+    M, N = size(X)
+
+    # Column synthesis
+    PCe = P.planθ[1]
+    PCo = P.planθ[2]
+
+    X[1] *= two(T)
+    A_mul_B_col_J!(Y, PCe, X, 1)
+    X[1] *= half(T)
+
+    for J = 2:4:N
+        A_mul_B_col_J!(Y, PCo, X, J)
+        A_mul_B_col_J!(Y, PCo, X, J+1)
+    end
+    for J = 4:4:N
+        X[1,J] *= two(T)
+        X[1,J+1] *= two(T)
+        A_mul_B_col_J!(Y, PCe, X, J)
+        A_mul_B_col_J!(Y, PCe, X, J+1)
+        X[1,J] *= half(T)
+        X[1,J+1] *= half(T)
+    end
+    scale!(half(T), Y)
+
+    # Row synthesis
+    scale!(inv(sqrt(π)), Y)
+    invsqrttwo = inv(sqrt(2))
+    @inbounds for i = 1:M Y[i] *= invsqrttwo end
+
+    temp = P.temp
+    planφ = P.planφ
+    C = P.C
+    for I = 1:M
+        copy_row_I!(temp, Y, I)
+        row_synthesis!(planφ, C, temp)
+        copy_row_I!(Y, temp, I)
+    end
+    Y
+end
+
+function Base.A_mul_B!{T}(Y::Matrix{T}, P::AnalysisPlan{T}, X::Matrix{T})
+    M, N = size(X)
+
+    # Row analysis
+    temp = P.temp
+    planφ = P.planφ
+    C = P.C
+    for I = 1:M
+        copy_row_I!(temp, X, I)
+        row_analysis!(planφ, C, temp)
+        copy_row_I!(Y, temp, I)
+    end
+
+    # Column analysis
+    PCe = P.planθ[1]
+    PCo = P.planθ[2]
+
+    A_mul_B_col_J!(Y, PCe, Y, 1)
+    Y[1] *= half(T)
+    for J = 2:4:N
+        A_mul_B_col_J!(Y, PCo, Y, J)
+        A_mul_B_col_J!(Y, PCo, Y, J+1)
+    end
+    for J = 4:4:N
+        A_mul_B_col_J!(Y, PCe, Y, J)
+        A_mul_B_col_J!(Y, PCe, Y, J+1)
+        Y[1,J] *= half(T)
+        Y[1,J+1] *= half(T)
+    end
+    scale!(sqrt(π)*inv(T(M)), Y)
+    sqrttwo = sqrt(2)
+    @inbounds for i = 1:M Y[i] *= sqrttwo end
+
+    Y
+end
+
+
+
+
+function row_analysis!{T}(P, C, vals::Vector{T})
+    n = length(vals)
+    cfs = scale!(two(T)/n,P*vals)
+    cfs[1] *= half(T)
+    if iseven(n)
+        cfs[n÷2+1] *= half(T)
+    end
+
+    negateeven!(reverseeven!(A_mul_B!(C, cfs)))
+end
+
+function row_synthesis!{T}(P, C, cfs::Vector{T})
+    n = length(cfs)
+    Ac_mul_B!(C, reverseeven!(negateeven!(cfs)))
+    if iseven(n)
+        cfs[n÷2+1] *= two(T)
+    end
+    cfs[1] *= two(T)
+    P*scale!(half(T), cfs)
+end
+
+function copy_row_I!(temp::Vector, Y::Matrix, I::Int)
+    M, N = size(Y)
+    @inbounds @simd for j = 1:N
+        temp[j] = Y[I+M*(j-1)]
+    end
+    temp
+end
+
+function copy_row_I!(Y::Matrix, temp::Vector, I::Int)
+    M, N = size(Y)
+    @inbounds @simd for j = 1:N
+        Y[I+M*(j-1)] = temp[j]
+    end
+    Y
+end
+
+
+function reverseeven!(x::Vector)
+    n = length(x)
+    if iseven(n)
+        @inbounds @simd for k=2:2:n÷2
+            x[k], x[n+2-k] = x[n+2-k], x[k]
+        end
+    else
+        @inbounds @simd for k=2:2:n÷2
+            x[k], x[n+1-k] = x[n+1-k], x[k]
+        end
+    end
+    x
+end
+
+function negateeven!(x::Vector)
+    @inbounds @simd for k = 2:2:length(x)
+        x[k] *= -1
+    end
+    x
+end
+
+import Base.FFTW: unsafe_execute!, fftwSingle, fftwDouble, fftwNumber
+import Base.FFTW: libfftw, libfftwf, PlanPtr, r2rFFTWPlan
+
+function A_mul_B_col_J!{T}(Y::Matrix{T}, P::r2rFFTWPlan{T}, X::Matrix{T}, J::Int)
+    unsafe_execute_col_J!(P, X, Y, J)
+    return Y
+end
+
+function unsafe_execute_col_J!{T<:fftwDouble}(plan::r2rFFTWPlan{T}, X::Matrix{T}, Y::Matrix{T}, J::Int)
+    M = size(X, 1)
+    ccall((:fftw_execute_r2r, libfftw), Void, (PlanPtr, Ptr{T}, Ptr{T}), plan, pointer(X, M*(J-1)+1), pointer(Y, M*(J-1)+1))
+end
+
+function unsafe_execute_col_J!{T<:fftwSingle}(plan::r2rFFTWPlan{T}, X::Matrix{T}, Y::Matrix{T}, J::Int)
+    M = size(X, 1)
+    ccall((:fftwf_execute_r2r, libfftwf), Void, (PlanPtr, Ptr{T}, Ptr{T}), plan, pointer(X, M*(J-1)+1), pointer(Y, M*(J-1)+1))
+end

test/sphericalharmonictests.jl

@@ -54,6 +54,12 @@ for θ in (0.123, 0.456), m in 2:2:n
     @test norm(S1-SA) < 1000eps()
 end
 
+println()
+println("Testing synthesis and analysis")
+println()
+
+include("synthesisanalysistests.jl")
+
 println()
 println("Testing API")
 println()

test/synthesisanalysistests.jl

@@ -0,0 +1,48 @@
+using FastTransforms, Base.Test
+
+import FastTransforms: normalizecolumns!, maxcolnorm
+
+# Starting with normalized spherical harmonic coefficients,
+
+n = 50
+F = sphrandn(Float64, n, n);
+normalizecolumns!(F);
+
+# we convert to bivariate Fourier series.
+
+G = sph2fourier(F);
+
+# At equispaced points in angle,
+
+θ = (0.5:n-0.5)*π/n
+φ = (0:2n-2)*2π/(2n-1)
+
+# the Fourier series evaluates to:
+
+SF = [sum(G[ℓ,1]*cos((ℓ-1)*θ)/sqrt(2π) for ℓ in 1:n) + sum(G[ℓ,2m]*sin(ℓ*θ)*sin(m*φ)/sqrt(π) for ℓ in 1:n, m in 1:2:n-1) + sum(G[ℓ,2m+1]*sin(ℓ*θ)*cos(m*φ)/sqrt(π) for ℓ in 1:n, m in 1:2:n-1) + sum(G[ℓ,2m]*cos((ℓ-1)*θ)*sin(m*φ)/sqrt(π) for ℓ in 1:n, m in 2:2:n-1) + sum(G[ℓ,2m+1]*cos((ℓ-1)*θ)*cos(m*φ)/sqrt(π) for ℓ in 1:n, m in 2:2:n-1) for θ in θ, φ in φ]
+
+# but that was slow, so we accelerate it via in-place FFTW technology:
+
+Ps = FastTransforms.plan_synthesis(G);
+
+Y = zero(G);
+
+A_mul_B!(Y, Ps, G)
+
+@test maxcolnorm(SF - Y) < 10000eps()
+
+# Retracing our steps, function values on the sphere are converted to Fourier coefficients:
+
+Pa = FastTransforms.plan_analysis(Y);
+
+Z = zero(Y);
+
+A_mul_B!(Z, Pa, Y)
+
+@test maxcolnorm(Z - G) < 10eps()
+
+# And Fourier coefficients are converted back to spherical harmonic coefficients:
+
+H = fourier2sph(Z)
+
+@test maxcolnorm(F - H) < 100eps()