Merge pull request #47 from andreasnoack/an/fix01

Fix the size one case and provide informative error message in the size zero case.

Author: dannys4

src/algos.jl

@@ -73,7 +73,7 @@ function fft_dft!(out::AbstractVector{T}, in::AbstractVector{T}, N::Int, start_o
         tmp += in[start_in + j*stride_in]
     end
     out[start_out] = tmp
-    
+
     wk = wkn = w = convert(T, cispi(direction_sign(d)*2/N))
     @inbounds for d in 1:N-1
         tmp = in[start_in]
@@ -98,7 +98,7 @@ function fft_dft!(out::AbstractVector{Complex{T}}, in::AbstractVector{T}, N::Int
     end
     out[start_out] = convert(Complex{T}, tmpBegin)
     iseven(N) && (out[start_out + stride_out*halfN] = convert(Complex{T}, tmpHalf))
-    
+
     @inbounds for d in 1:halfN
         tmp = in[start_in]
         @inbounds for k in 1:N-1
@@ -200,7 +200,7 @@ function fft_pow4!(out::AbstractVector{T}, in::AbstractVector{U}, N::Int, start_
 
     w1 = convert(T, cispi(direction_sign(d)*2/N))
     wj = one(T)
-    
+
     w1 = convert(T, cispi(ds*2/N))
     w2 = convert(T, cispi(ds*4/N))
     w3 = convert(T, cispi(ds*6/N))
@@ -212,7 +212,7 @@ function fft_pow4!(out::AbstractVector{T}, in::AbstractVector{U}, N::Int, start_
         @muladd k2 = start_out + (k+2*m)*stride_out
         @muladd k3 = start_out + (k+3*m)*stride_out
         y_k0, y_k1, y_k2, y_k3 = out[k0], out[k1], out[k2], out[k3]
-        @muladd out[k0] = (y_k0 + y_k2*wk2) + (y_k1*wk1 + y_k3*wk2)
+        @muladd out[k0] = (y_k0 + y_k2*wk2) + (y_k1*wk1 + y_k3*wk3)
         @muladd out[k1] = (y_k0 - y_k2*wk2) + (y_k1*wk1 - y_k3*wk3) * plusi
         @muladd out[k2] = (y_k0 + y_k2*wk2) - (y_k1*wk1 + y_k3*wk3)
         @muladd out[k3] = (y_k0 - y_k2*wk2) + (y_k1*wk1 - y_k3*wk3) * minusi

src/callgraph.jl

@@ -55,18 +55,6 @@ end
 # Get the node in the graph at index i
 Base.getindex(g::CallGraph{T}, i::Int) where {T} = g.nodes[i]
 
-"""
-$(TYPEDSIGNATURES)
-Check if `N` is a power of `base`
-
-"""
-function _ispow(N, base)
-    while N % base == 0
-        N = N/base
-    end
-    return N == 1
-end
-
 """
 $(TYPEDSIGNATURES)
 Check if `N` is a power of 2 or 4
@@ -93,6 +81,9 @@ Recursively instantiate a set of `CallGraphNode`s
 
 """
 function CallGraphNode!(nodes::Vector{CallGraphNode}, N::Int, workspace::Vector{Vector{T}}, s_in::Int, s_out::Int)::Int where {T}
+    if N == 0
+        throw(DimensionMismatch("array has to be non-empty"))
+    end
     if iseven(N)
         pow = _ispow24(N)
         if !isnothing(pow)
@@ -102,15 +93,15 @@ function CallGraphNode!(nodes::Vector{CallGraphNode}, N::Int, workspace::Vector{
         end
     end
     if N % 3 == 0
-        if _ispow(N, 3)
+        if nextpow(3, N) == N
             push!(workspace, T[])
             push!(nodes, CallGraphNode(0, 0, Pow3FFT(), N, s_in, s_out))
             return 1
         end
     end
-    if isprime(N)
+    if N == 1 || isprime(N)
         push!(workspace, T[])
-        push!(nodes, CallGraphNode(0,0, DFT(),N, s_in, s_out))
+        push!(nodes, CallGraphNode(0, 0, DFT(), N, s_in, s_out))
         return 1
     end
     Ns = [first(x) for x in collect(factor(N)) for _ in 1:last(x)]

test/onedim/complex_backward.jl

@@ -8,4 +8,13 @@ test_nums = [8, 11, 15, 16, 27, 100]
         y_ref[1] = N
         @test y ≈ y_ref atol=1e-12
     end
-end
+end
+
+@testset verbose = true "against naive implementation. Size: $n" for n in 1:64
+    x = complex.(randn(n), randn(n))
+    @test naive_1d_fourier_transform(x, FFTA.FFT_BACKWARD) ≈ bfft(x)
+end
+
+@testset "error messages" begin
+    @test_throws DimensionMismatch bfft(zeros(0))
+end

test/onedim/complex_forward.jl

@@ -8,4 +8,13 @@ test_nums = [8, 11, 15, 16, 27, 100]
         y_ref[1] = N
         @test y ≈ y_ref atol=1e-12
     end
-end
+end
+
+@testset verbose = true "against naive implementation. Size: $n" for n in 1:64
+    x = complex.(randn(n), randn(n))
+    @test naive_1d_fourier_transform(x, FFTA.FFT_FORWARD) ≈ fft(x)
+end
+
+@testset "error messages" begin
+    @test_throws DimensionMismatch fft(zeros(0))
+end

test/onedim/real_backward.jl

@@ -11,4 +11,20 @@ test_nums = [8, 11, 15, 16, 27, 100]
         end
         @test y_ref ≈ y atol=1e-12
     end
-end
+end
+
+@testset verbose = true "against naive implementation. Size: $n" for n in 1:64
+    x = complex.(randn(n ÷ 2 + 1), randn(n ÷ 2 + 1))
+    x[begin] = real(x[begin])
+    if iseven(n)
+        x[end] = real(x[end])
+        xe = [x; conj.(reverse(x[begin + 1:end - 1]))]
+    else
+        xe = [x; conj.(reverse(x[begin + 1:end]))]
+    end
+    @test naive_1d_fourier_transform(xe, FFTA.FFT_BACKWARD) ≈ brfft(x, n)
+end
+
+@testset "error messages" begin
+    @test_throws DimensionMismatch brfft(zeros(ComplexF64, 0), 0)
+end

test/onedim/real_forward.jl

@@ -1,11 +1,20 @@
 using FFTA, Test
 test_nums = [8, 11, 15, 16, 27, 100]
-@testset verbose = true " forward" begin 
+@testset verbose = true " forward" begin
     for N in test_nums
         x = ones(Float64, N)
         y = rfft(x)
         y_ref = 0*y
         y_ref[1] = N
         @test y ≈ y_ref atol=1e-12
     end
+end
+
+@testset verbose = true "against naive implementation. Size: $n" for n in 1:64
+    x = randn(n)
+    @test naive_1d_fourier_transform(x, FFTA.FFT_FORWARD)[1:(n ÷ 2 + 1)] ≈ rfft(x)
+end
+
+@testset "error messages" begin
+    @test_throws DimensionMismatch rfft(zeros(0))
 end

test/runtests.jl

@@ -1,4 +1,4 @@
-using Test, Random
+using Test, Random, FFTA
 
 function padnum(m,x)
     digs = floor(Int, log10(m))
@@ -8,6 +8,42 @@ function padnum(m,x)
     String(v)
 end
 
+function naive_1d_fourier_transform(x::Vector, d::FFTA.Direction)
+    n = length(x)
+    y = zeros(Complex{Float64}, n)
+
+    for u in 0:(n - 1)
+        s = 0.0 + 0.0im
+        for v in 0:(n - 1)
+            a = FFTA.direction_sign(d) * 2π * u * v / n
+            s += x[v + 1] * exp(im * a)
+        end
+        y[u + 1] = s
+    end
+
+    return y
+end
+
+function naive_2d_fourier_transform(X::Matrix, d::FFTA.Direction)
+    rows, cols = size(X)
+    Y = zeros(Complex{Float64}, rows, cols)
+
+    for u in 0:(rows - 1)
+        for v in 0:(cols - 1)
+            s = 0.0 + 0.0im
+            for x in 0:(rows - 1)
+                for y in 0:(cols - 1)
+                    a = FFTA.direction_sign(d) * 2π * (u * x / rows + v * y / cols)
+                    s += X[x + 1, y + 1] * exp(im * a)
+                end
+            end
+            Y[u + 1, v + 1] = s
+        end
+    end
+
+    return Y
+end
+
 Random.seed!(1)
 @testset verbose = true "FFTA" begin
     @testset verbose = true "1D" begin

test/twodim/complex_backward.jl

@@ -8,4 +8,16 @@ test_nums = [8, 11, 15, 16, 27, 100]
         y_ref[1] = length(x)
         @test y ≈ y_ref
     end
-end
+end
+
+@testset verbose = true "against naive implementation" for n in 1:64
+    @testset "size: ($m, $n)" for m in n:(n + 1)
+        X = complex.(randn(m, n), randn(m, n))
+        Y = similar(X)
+        @test naive_2d_fourier_transform(X, FFTA.FFT_BACKWARD) ≈ bfft(X)
+    end
+end
+
+@testset "error messages" begin
+    @test_throws DimensionMismatch bfft(zeros(ComplexF64, 0, 0))
+end

test/twodim/complex_forward.jl

@@ -8,4 +8,16 @@ test_nums = [8, 11, 15, 16, 27, 100]
         y_ref[1] = length(x)
         @test y ≈ y_ref
     end
-end
+end
+
+@testset verbose = true "against naive implementation" for n in 1:64
+    @testset "size: ($m, $n)" for m in n:(n + 1)
+        X = complex.(randn(m, n), randn(m, n))
+        Y = similar(X)
+        @test naive_2d_fourier_transform(X, FFTA.FFT_FORWARD) ≈ fft(X)
+    end
+end
+
+@testset "error messages" begin
+    @test_throws DimensionMismatch fft(zeros(0, 0))
+end