Bifrost FFT != numpy

[telegraphic]

This fails:

""" Test of FFT bifrost kernel """

from bifrost.libbifrost import _bf
from bifrost.GPUArray import GPUArray
import bifrost as bf
import numpy as np
from bifrost.fft import fft

np.set_printoptions(precision=2)

n_iter = 1    # Number of iterations to run

for dx in (1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 20000000):

    print "\nFFT SIZE: %i" % dx
    print "-------------"

    # Setup some input / output arrays
    input_data   = GPUArray(shape=(dx), dtype=np.complex64)
    output_data  = GPUArray(shape=(dx), dtype=np.complex64)

    tmp_arr = np.arange(dx, dtype=np.complex64)
    input_data.set(tmp_arr)

    # Run the accumulator
    print "Comparing FFT"
    for ii in range(n_iter):
        print "%i of %i"% (ii+1, n_iter)
        err = fft(input_data.as_BFarray(), output_data.as_BFarray())
        dd = output_data.get()
        zz = np.fft.fft(tmp_arr)
        try:
            assert np.allclose(dd, zz, rtol=0.01)
        except:
            deltas = np.abs(dd - zz)
            print "ERROR: data does not match"
            print "Error code from FFT: %s" % err
            print "Max delta: %2.2fs @ index: %i" % (np.max(deltas), np.argmax(deltas))

            q = np.isclose(dd, zz, rtol=0.01)
            fail_idx = np.where(q==False)[0]
            print "Number elements failing test: %s" % len(fail_idx)
            print "Failing indexes:"
            print fail_idx
            print "Bifrost: ", dd[fail_idx]
            print "Numpy  : ", zz[fail_idx]

with this output:

FFT SIZE: 10000
-------------
Comparing FFT
1 of 1

FFT SIZE: 100000
-------------
Comparing FFT
1 of 1

FFT SIZE: 1000000
-------------
Comparing FFT
1 of 1
ERROR: data does not match
Error code from FFT: 0
Max delta: 34462.19s @ index: 999999
Number elements failing test: 10
Failing indexes:
[200000 240000 376000 495999 499200 500800 559999 624000 760000 800000]
Bifrost:  [-515097.72+696916.5j  -500086.81+544187.62j -506977.94+204702.j
 -493507.59  +3326.22j -499968.16 +13998.7j  -499967.91 -13996.47j
 -494384.22 -92309.42j -506977.94-204701.86j -500086.78-544187.56j
 -515097.72-696916.5j ]
Numpy  :  [-500000.+688190.96j -500000.+532445.92j -500000.+205268.86j
 -500000.  +6285.09j -500000.  +1256.64j -500000.  -1256.64j
 -500000. -95378.47j -500000.-205268.86j -500000.-532445.92j
 -500000.-688190.96j]

FFT SIZE: 10000000
-------------
Comparing FFT
1 of 1
ERROR: data does not match
Error code from FFT: 0
Max delta: 5092706.40s @ index: 9999999
Number elements failing test: 1293
Failing indexes:
[ 240000  319999  400000 ..., 9376000 9600000 9760000]
Bifrost:  [-4842115.0+67201176.j -4701928.5+50007032.j -5212520.5+39969704.j ...,
 -4857975.5-25402098.j -5212520.0-39969704.j -4842115.0-67201184.j]
Numpy  :  [-5000000.+66188848.26j -5000000.+49568411.06j -5000000.+39579075.44j ...,
 -5000000.-25178034.07j -5000000.-39579075.44j -5000000.-66188848.26j]

Although TBH it looks like numpy might be wrong?

[benbarsdell]

My guess is a precision issue due to pathological input data.

If I replace the input data with this instead: tmp_arr = np.random.random(size=dx*2).astype(np.float32).view(np.complex64) then the tests all pass.

Also worth noting that Numpy only supports double-precision FFTs (the output is always complex128). As an alternative, scipy.fftpack.fft supports single-precision.

[telegraphic]

Interesting... OK, closing this for now