add test for uniform points/frequencies

- change `cg` to `cg_for_inufft` in case the package is imported. (It’s
less likely that this name would be used elsewhere.)

- add `cg` vectors to the `iNUFFTPlan`

Author: MikaelSlevinsky

src/inufft.jl

@@ -6,6 +6,9 @@ For best performance, choose the right number of threads by `FFTW.set_num_thread
 immutable iNUFFTPlan{N,T,S,PT} <: Base.DFT.Plan{T}
     pt::PT
     TP::Toeplitz{T}
+    r::Vector{T}
+    p::Vector{T}
+    Ap::Vector{T}
     ϵ::S
 end
 
@@ -22,8 +25,11 @@ function plan_inufft1{T<:AbstractFloat}(ω::AbstractVector{T}, ϵ::T)
     r[1] = avg
     c[1] = avg
     TP = Toeplitz(c, r)
+    r = zero(c)
+    p = zero(c)
+    Ap = zero(c)
 
-    iNUFFTPlan{1, eltype(TP), typeof(ϵ), typeof(pt)}(pt, TP, ϵ)
+    iNUFFTPlan{1, eltype(TP), typeof(ϵ), typeof(pt)}(pt, TP, r, p, Ap, ϵ)
 end
 
 doc"""
@@ -38,23 +44,26 @@ function plan_inufft2{T<:AbstractFloat}(x::AbstractVector{T}, ϵ::T)
     r[1] = avg
     c[1] = avg
     TP = Toeplitz(c, r)
+    r = zero(c)
+    p = zero(c)
+    Ap = zero(c)
 
-    iNUFFTPlan{2, eltype(TP), typeof(ϵ), typeof(pt)}(pt, TP, ϵ)
+    iNUFFTPlan{2, eltype(TP), typeof(ϵ), typeof(pt)}(pt, TP, r, p, Ap, ϵ)
 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ϵ)
+    pt, TP, r, p, Ap, ϵ = P.pt, P.TP, P.r, P.p, P.Ap, P.ϵ
+    cg_for_inufft(TP, c, f, r, p, Ap, 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ϵ)
+    pt, TP, r, p, Ap, ϵ = P.pt, P.TP, P.r, P.p, P.Ap, P.ϵ
+    cg_for_inufft(TP, c, conj!(pt*conj!(f)), r, p, Ap, 50, 100ϵ)
     conj!(f)
     c
 end
@@ -69,16 +78,16 @@ 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)
+function cg_for_inufft{T}(A::ToeplitzMatrices.AbstractToeplitz{T}, x::AbstractVector{T}, b::AbstractVector{T}, r::AbstractVector{T}, p::AbstractVector{T}, Ap::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)
+    fill!(r, zero(T))
+    fill!(p, zero(T))
+    fill!(Ap, zero(T))
     # r = b - A*x
     copy!(r, b)
     A_mul_B!(-one(T), A, x, one(T), r)

src/nufft.jl

@@ -301,7 +301,7 @@ F_{i,j} = \sum_{k=0}^{M-1}\sum_{\ell=0}^{N-1} C_{k,\ell} e^{-2\pi{\rm i} (x_i k
 nufft2{T<:AbstractFloat}(C::Matrix, x::AbstractVector{T}, y::AbstractVector{T}, ϵ::T) = plan_nufft2(x, y, ϵ)*C
 
 
-FindK{T<:AbstractFloat}(γ::T, ϵ::T) = Int(ceil(5*γ*exp(lambertw(log(10/ϵ)/γ/7))))
+FindK{T<:AbstractFloat}(γ::T, ϵ::T) = γ < ϵ ? 1 : Int(ceil(5*γ*exp(lambertw(log(10/ϵ)/γ/7))))
 AssignClosestEquispacedGridpoint{T<:AbstractFloat}(x::AbstractVector{T})::AbstractVector{T} = round.([Int], size(x, 1)*x)
 AssignClosestEquispacedFFTpoint{T<:AbstractFloat}(x::AbstractVector{T})::Array{Int,1} = mod.(round.([Int], size(x, 1)*x), size(x, 1)) + 1
 PerturbationParameter{T<:AbstractFloat}(x::AbstractVector{T}, s_vec::AbstractVector{T})::AbstractFloat = norm(size(x, 1)*x - s_vec, Inf)

test/nuffttests.jl

@@ -68,6 +68,18 @@ using FastTransforms, Base.Test
         @test norm(exact - fast, Inf) < err_bnd
     end
 
+    # Check that if points/frequencies are indeed uniform, then it's equal to the fft.
+
+    n = 1000
+    c = complex(rand(n))
+    ω = collect(0.0:n-1)
+    x = ω/n
+    fftc = fft(c)
+    @test norm(nufft1(c, ω, eps()) - fftc) == 0
+    @test norm(nufft2(c, x, eps()) - fftc) == 0
+    @test norm(nufft3(c, x, ω, eps()) - fftc) == 0
+
+
     function nudft1{T<:AbstractFloat}(C::Matrix{Complex{T}}, ω1::AbstractVector{T}, ω2::AbstractVector{T})
         # Nonuniform discrete Fourier transform of type I-I