Merge pull request #22 from MikaelSlevinsky/add-nufft

Add nonuniform fast Fourier transforms

Author: MikaelSlevinsky

README.md

@@ -75,6 +75,48 @@ is valid for the half-open square `(α,β) ∈ (-1/2,1/2]^2`. Therefore, the fas
 when the parameters are inside. If the parameters `(α,β)` are not exceptionally beyond the square,
 then increment/decrement operators are used with linear complexity (and linear conditioning) in the degree.
 
+## Nonuniform fast Fourier transforms
+
+The NUFFTs are implemented thanks to [Alex Townsend](https://github.com/ajt60gaibb):
+ - `nufft1` assumes uniform samples and noninteger frequencies;
+ - `nufft2` assumes nonuniform samples and integer frequencies;
+ - `nufft3 ( = nufft)` assumes nonuniform samples and noninteger frequencies;
+ - `inufft1` inverts an `nufft1`; and,
+ - `inufft2` inverts an `nufft2`.
+
+Here is an example:
+```julia
+julia> n = 10^4;
+
+julia> c = complex(rand(n));
+
+julia> ω = collect(0:n-1) + rand(n);
+
+julia> nufft1(c, ω, eps());
+
+julia> p1 = plan_nufft1(ω, eps());
+
+julia> @time p1*c;
+  0.002383 seconds (6 allocations: 156.484 KiB)
+
+julia> x = (collect(0:n-1) + 3rand(n))/n;
+
+julia> nufft2(c, x, eps());
+
+julia> p2 = plan_nufft2(x, eps());
+
+julia> @time p2*c;
+  0.001478 seconds (6 allocations: 156.484 KiB)
+
+julia> nufft3(c, x, ω, eps());
+
+julia> p3 = plan_nufft3(x, ω, eps());
+
+julia> @time p3*c;
+  0.058999 seconds (6 allocations: 156.484 KiB)
+
+```
+
 ## The Padua Transform
 
 The Padua transform and its inverse are implemented thanks to [Michael Clarke](https://github.com/MikeAClarke). These are optimized methods designed for computing the bivariate Chebyshev coefficients by interpolating a bivariate function at the Padua points on `[-1,1]^2`.
@@ -87,7 +129,7 @@ julia> N = div((n+1)*(n+2),2);
 julia> v = rand(N); # The length of v is the number of Padua points
 
 julia> @time norm(ipaduatransform(paduatransform(v))-v)
-0.006571 seconds (846 allocations: 1.746 MiB)
+  0.006571 seconds (846 allocations: 1.746 MiB)
 3.123637691861415e-14
 
 ```
@@ -126,10 +168,12 @@ As with other fast transforms, `plan_sph2fourier` saves effort by caching the pr
 
    [2]  N. Hale and A. Townsend. <a href="http://dx.doi.org/10.1137/130932223">A fast, simple, and stable Chebyshev—Legendre transform using and asymptotic formula</a>, *SIAM J. Sci. Comput.*, **36**:A148—A167, 2014.
 
-   [3] J. Keiner. <a href="http://dx.doi.org/10.1137/070703065">Computing with expansions in Gegenbauer polynomials</a>, *SIAM J. Sci. Comput.*, **31**:2151—2171, 2009.
+   [3]  J. Keiner. <a href="http://dx.doi.org/10.1137/070703065">Computing with expansions in Gegenbauer polynomials</a>, *SIAM J. Sci. Comput.*, **31**:2151—2171, 2009.
+
+   [4]  D. Ruiz—Antolín and A. Townsend. <a href="https://arxiv.org/abs/1701.04492">A nonuniform fast Fourier transform based on low rank approximation</a>, arXiv:1701.04492, 2017.
 
-   [4]  R. M. Slevinsky. <a href="https://doi.org/10.1093/imanum/drw070">On the use of Hahn's asymptotic formula and stabilized recurrence for a fast, simple, and stable Chebyshev—Jacobi transform</a>, in press at *IMA J. Numer. Anal.*, 2017.
+   [5]  R. M. Slevinsky. <a href="https://doi.org/10.1093/imanum/drw070">On the use of Hahn's asymptotic formula and stabilized recurrence for a fast, simple, and stable Chebyshev—Jacobi transform</a>, in press at *IMA J. Numer. Anal.*, 2017.
 
-   [5]  R. M. Slevinsky. <a href="https://arxiv.org/abs/1705.05448">Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series</a>, arXiv:1705.05448, 2017.
+   [6]  R. M. Slevinsky. <a href="https://arxiv.org/abs/1705.05448">Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series</a>, arXiv:1705.05448, 2017.
 
-   [6]  A. Townsend, M. Webb, and S. Olver. <a href="https://doi.org/10.1090/mcom/3277">Fast polynomial transforms based on Toeplitz and Hankel matrices</a>, in press at *Math. Comp.*, 2017.
+   [7]  A. Townsend, M. Webb, and S. Olver. <a href="https://doi.org/10.1090/mcom/3277">Fast polynomial transforms based on Toeplitz and Hankel matrices</a>, in press at *Math. Comp.*, 2017.

docs/src/index.md

@@ -52,6 +52,46 @@ plan_cjt
 plan_icjt
 ```
 
+```@docs
+nufft1
+```
+
+```@docs
+nufft2
+```
+
+```@docs
+nufft3
+```
+
+```@docs
+inufft1
+```
+
+```@docs
+inufft2
+```
+
+```@docs
+plan_nufft1
+```
+
+```@docs
+plan_nufft2
+```
+
+```@docs
+plan_nufft3
+```
+
+```@docs
+plan_inufft1
+```
+
+```@docs
+plan_inufft2
+```
+
 ```@docs
 paduatransform
 ```
@@ -112,6 +152,10 @@ FastTransforms.δ
 FastTransforms.Λ
 ```
 
+```@docs
+FastTransforms.lambertw
+```
+
 ```@docs
 FastTransforms.pochhammer
 ```

src/FastTransforms.jl

@@ -16,7 +16,12 @@ export normleg2cheb, cheb2normleg, normleg12cheb2, cheb22normleg1
 export plan_leg2cheb, plan_cheb2leg
 export plan_normleg2cheb, plan_cheb2normleg
 export plan_normleg12cheb2, plan_cheb22normleg1
+
 export gaunt
+
+export nufft, nufft1, nufft2, nufft3, inufft1, inufft2
+export plan_nufft, plan_nufft1, plan_nufft2, plan_nufft3, plan_inufft1, plan_inufft2
+
 export paduatransform, ipaduatransform, paduatransform!, ipaduatransform!, paduapoints
 export plan_paduatransform!, plan_ipaduatransform!
 
@@ -49,6 +54,8 @@ include("cheb2jac.jl")
 include("ChebyshevUltrasphericalPlan.jl")
 include("ultra2cheb.jl")
 include("cheb2ultra.jl")
+include("nufft.jl")
+include("inufft.jl")
 
 include("cjt.jl")
 

src/inufft.jl

@@ -0,0 +1,104 @@
+doc"""
+Pre-computes an inverse nonuniform fast Fourier transform of type `N`.
+
+For best performance, choose the right number of threads by `FFTW.set_num_threads(4)`, for example.
+"""
+immutable iNUFFTPlan{N,T,S,PT} <: Base.DFT.Plan{T}
+    pt::PT
+    TP::Toeplitz{T}
+    ϵ::S
+end
+
+doc"""
+Pre-computes an inverse nonuniform fast Fourier transform of type I.
+"""
+function plan_inufft1{T<:AbstractFloat}(ω::AbstractVector{T}, ϵ::T)
+    N = length(ω)
+    p = plan_nufft1(ω, ϵ)
+    pt = plan_nufft2(ω/N, ϵ)
+    c = p*ones(Complex{T}, N)
+    r = conj(c)
+    avg = (r[1]+c[1])/2
+    r[1] = avg
+    c[1] = avg
+    TP = Toeplitz(c, r)
+
+    iNUFFTPlan{1, eltype(TP), typeof(ϵ), typeof(pt)}(pt, TP, ϵ)
+end
+
+doc"""
+Pre-computes an inverse nonuniform fast Fourier transform of type II.
+"""
+function plan_inufft2{T<:AbstractFloat}(x::AbstractVector{T}, ϵ::T)
+    N = length(x)
+    pt = plan_nufft1(N*x, ϵ)
+    r = pt*ones(Complex{T}, N)
+    c = conj(r)
+    avg = (r[1]+c[1])/2
+    r[1] = avg
+    c[1] = avg
+    TP = Toeplitz(c, r)
+
+    iNUFFTPlan{2, eltype(TP), typeof(ϵ), typeof(pt)}(pt, TP, ϵ)
+end
+
+function (*){N,T,V}(p::iNUFFTPlan{N,T}, x::AbstractVector{V})
+    A_mul_B!(zeros(promote_type(T,V), length(x)), p, x)
+end
+
+function Base.A_mul_B!{T}(c::AbstractVector{T}, P::iNUFFTPlan{1,T}, f::AbstractVector{T})
+    pt, TP, ϵ = P.pt, P.TP, P.ϵ
+    cg(TP, c, f, 50, 100ϵ)
+    conj!(A_mul_B!(c, pt, conj!(c)))
+end
+
+function Base.A_mul_B!{T}(c::AbstractVector{T}, P::iNUFFTPlan{2,T}, f::AbstractVector{T})
+    pt, TP, ϵ = P.pt, P.TP, P.ϵ
+    cg(TP, c, conj!(pt*conj!(f)), 50, 100ϵ)
+    conj!(f)
+    c
+end
+
+doc"""
+Computes an inverse nonuniform fast Fourier transform of type I.
+"""
+inufft1{T<:AbstractFloat}(c::AbstractVector, ω::AbstractVector{T}, ϵ::T) = plan_inufft1(ω, ϵ)*c
+
+doc"""
+Computes an inverse nonuniform fast Fourier transform of type II.
+"""
+inufft2{T<:AbstractFloat}(c::AbstractVector, x::AbstractVector{T}, ϵ::T) = plan_inufft2(x, ϵ)*c
+
+function cg{T}(A::ToeplitzMatrices.AbstractToeplitz{T}, x::AbstractVector{T}, b::AbstractVector{T}, max_it::Integer, tol::Real)
+	n = length(b)
+	n1, n2 = size(A)
+	n == n1 == n2 || throw(DimensionMismatch(""))
+    nrmb = norm(b)
+    if nrmb == 0 nrmb = one(typeof(nrmb)) end
+	copy!(x, b)
+    r = zero(x)
+    p = zero(x)
+    Ap = zero(x)
+    # r = b - A*x
+    copy!(r, b)
+    A_mul_B!(-one(T), A, x, one(T), r)
+	copy!(p, r)
+	nrm2 = r⋅r
+    for k = 1:max_it
+        # Ap = A*p
+        A_mul_B!(one(T), A, p, zero(T), Ap)
+		α = nrm2/(p⋅Ap)
+        @inbounds @simd for l = 1:n
+            x[l] += α*p[l]
+            r[l] -= α*Ap[l]
+        end
+		nrm2new = r⋅r
+        cst = nrm2new/nrm2
+        @inbounds @simd for l = 1:n
+            p[l] = muladd(cst, p[l], r[l])
+        end
+		nrm2 = nrm2new
+        if sqrt(abs(nrm2)) ≤ tol*nrmb break end
+	end
+    return x
+end