[SciPy-dev] OS X fft problem is serious

[Tom Loredo]

Hi folks,

I hate to whine, especially when I don't know how to go about
solving the problem I'm whining about, but I don't think
it's good "advertising" to base the first official release of
new SciPy on a revision with a badly broken fft exposed in
the top-level namespace on widely-used platform.

Chris Fonnesbeck has verified that the problem exists also on
his Tiger build.  This is not a quirk of my Panther installation.
Something is really broken.

The ipython session logged below puts it starkly.  There is
a simple sine-like input array:  [0 1 0 -1 0].  fft and then
ifft with numpy and scipy.  numpy gives the array back (line 4).
scipy gives nonsense (line 7).

There is a clue to what is going on.  The ffts, though not
exactly equal, are essentially equal before the ifft (line 6).
But *after* taking the ifft, scipy changes the input fft
in-place (compare line 10 and line 9).  I believe the
default behavior is no overwrite.  Perhaps this points toward
a solution.  Perhaps the problem is not with SciPy but with
FFTW2; either way, I think we need to figure it out, preferably
before announcing the new release.

I wish I could help more with this, but I'll be at a workshop
for the next ~2 weeks.  I'll be glad to help when I'm free again.

-Tom

In [1]: import numpy, scipy
Overwriting fft=<function fft at 0x74d4f0> from scipy.fftpack.basic (was <function fft at 0x73c4f0> from 
numpy.dft.fftpack)
Overwriting ifft=<function ifft at 0x74d530> from scipy.fftpack.basic (was <function inverse_fft at 0x73c530> from 
numpy.dft.fftpack)

In [2]: a=numpy.array([0.,1.,0.,-1.,0.])

In [3]: nfa=numpy.fft(a)

In [4]: print numpy.ifft(nfa).real
[  0.00000000e+00   1.00000000e+00   9.76996262e-16  -1.00000000e+00
  -9.76996262e-16]

In [5]: sfa=scipy.fft(a)

In [6]: numpy.allclose(nfa,sfa)
Out[6]: True

In [7]: print scipy.ifft(sfa).real
[ 0.   0.5 -0.5 -0.5  0.5]

In [8]: numpy.allclose(nfa,sfa)
Out[8]: False

In [9]: print nfa
[ 0.        +0.j          1.11803399-1.53884177j -1.11803399+0.36327126j
 -1.11803399-0.36327126j  1.11803399+1.53884177j]

In [10]: print sfa
[ 0. +0.j  0.5+0.j -0.5-0.j -0.5+0.j  0.5-0.j]

In [11]: numpy.__version__
Out[11]: '0.9.3.1863'

In [12]: scipy.__version__
Out[12]: '0.4.4.1544'




-------------------------------------------------
This mail sent through IMP: http://horde.org/imp/