[Scipy-svn] r6570 - in branches/0.8.x: doc/release scipy/fftpack scipy/fftpack/tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sun Jun 27 08:53:15 EDT 2010
Author: ptvirtan
Date: 2010-06-27 07:53:15 -0500 (Sun, 27 Jun 2010)
New Revision: 6570
Modified:
branches/0.8.x/doc/release/0.8.0-notes.rst
branches/0.8.x/scipy/fftpack/basic.py
branches/0.8.x/scipy/fftpack/tests/test_basic.py
Log:
BUG: fftpack: use single-precision FFT implementation only for "easy" sizes
For large "difficult" sizes, the current single precision implementation
apparently falls back to a non-optimal algorithm, and can produce large
rounding errors. (bug #1212)
This commit provides a work-around, by enabling the single-precision
algorithms only for sizes that are composite numbers of 2, 3, and 5, for
which FFTPACK has explicit support. This is probably a rather
conservative approach.
Modified: branches/0.8.x/doc/release/0.8.0-notes.rst
===================================================================
--- branches/0.8.x/doc/release/0.8.0-notes.rst 2010-06-27 04:20:07 UTC (rev 6569)
+++ branches/0.8.x/doc/release/0.8.0-notes.rst 2010-06-27 12:53:15 UTC (rev 6570)
@@ -83,8 +83,7 @@
---------------------------
New realtransforms have been added, namely dct and idct for Discrete Cosine
-Transform; type I, II and III are available, for both single and double
-precision.
+Transform; type I, II and III are available.
Single precision support for fft functions (scipy.fftpack)
----------------------------------------------------------
@@ -92,6 +91,10 @@
fft functions can now handle single precision inputs as well: fft(x) will
return a single precision array if x is single precision.
+At the moment, for FFT sizes that are not composites of 2, 3, and 5, the
+transform is computed internally in double precision to avoid rounding error in
+FFTPACK.
+
Correlation functions now implement the usual definition (scipy.signal)
-----------------------------------------------------------------------
Modified: branches/0.8.x/scipy/fftpack/basic.py
===================================================================
--- branches/0.8.x/scipy/fftpack/basic.py 2010-06-27 04:20:07 UTC (rev 6569)
+++ branches/0.8.x/scipy/fftpack/basic.py 2010-06-27 12:53:15 UTC (rev 6570)
@@ -22,22 +22,73 @@
def istype(arr, typeclass):
return issubclass(arr.dtype.type, typeclass)
+# XXX: single precision FFTs partially disabled due to accuracy issues
+# for large prime-sized inputs.
+#
+# See http://permalink.gmane.org/gmane.comp.python.scientific.devel/13834
+# ("fftpack test failures for 0.8.0b1", Ralf Gommers, 17 Jun 2010,
+# @ scipy-dev)
+#
+# These should be re-enabled once the problems are resolved
+
+def _is_safe_size(n):
+ """
+ Is the size of FFT such that FFTPACK can handle it in single precision
+ with sufficient accuracy?
+
+ Composite numbers of 2, 3, and 5 are accepted, as FFTPACK has those
+ """
+ n = int(n)
+ for c in (2, 3, 5):
+ while n % c == 0:
+ n /= c
+ return (n <= 1)
+
+def _fake_crfft(x, n, *a, **kw):
+ if _is_safe_size(n):
+ return _fftpack.crfft(x, n, *a, **kw)
+ else:
+ return _fftpack.zrfft(x, n, *a, **kw).astype(numpy.complex64)
+
+def _fake_cfft(x, n, *a, **kw):
+ if _is_safe_size(n):
+ return _fftpack.cfft(x, n, *a, **kw)
+ else:
+ return _fftpack.zfft(x, n, *a, **kw).astype(numpy.complex64)
+
+def _fake_rfft(x, n, *a, **kw):
+ if _is_safe_size(n):
+ return _fftpack.rfft(x, n, *a, **kw)
+ else:
+ return _fftpack.drfft(x, n, *a, **kw).astype(numpy.float32)
+
+def _fake_cfftnd(x, shape, *a, **kw):
+ if numpy.all(map(_is_safe_size, shape)):
+ return _fftpack.cfftnd(x, shape, *a, **kw)
+ else:
+ return _fftpack.zfftnd(x, shape, *a, **kw).astype(numpy.complex64)
+
_DTYPE_TO_FFT = {
- numpy.dtype(numpy.float32): _fftpack.crfft,
+# numpy.dtype(numpy.float32): _fftpack.crfft,
+ numpy.dtype(numpy.float32): _fake_crfft,
numpy.dtype(numpy.float64): _fftpack.zrfft,
- numpy.dtype(numpy.complex64): _fftpack.cfft,
+# numpy.dtype(numpy.complex64): _fftpack.cfft,
+ numpy.dtype(numpy.complex64): _fake_cfft,
numpy.dtype(numpy.complex128): _fftpack.zfft,
}
_DTYPE_TO_RFFT = {
- numpy.dtype(numpy.float32): _fftpack.rfft,
+# numpy.dtype(numpy.float32): _fftpack.rfft,
+ numpy.dtype(numpy.float32): _fake_rfft,
numpy.dtype(numpy.float64): _fftpack.drfft,
}
_DTYPE_TO_FFTN = {
- numpy.dtype(numpy.complex64): _fftpack.cfftnd,
+# numpy.dtype(numpy.complex64): _fftpack.cfftnd,
+ numpy.dtype(numpy.complex64): _fake_cfftnd,
numpy.dtype(numpy.complex128): _fftpack.zfftnd,
- numpy.dtype(numpy.float32): _fftpack.cfftnd,
+# numpy.dtype(numpy.float32): _fftpack.cfftnd,
+ numpy.dtype(numpy.float32): _fake_cfftnd,
numpy.dtype(numpy.float64): _fftpack.zfftnd,
}
Modified: branches/0.8.x/scipy/fftpack/tests/test_basic.py
===================================================================
--- branches/0.8.x/scipy/fftpack/tests/test_basic.py 2010-06-27 04:20:07 UTC (rev 6569)
+++ branches/0.8.x/scipy/fftpack/tests/test_basic.py 2010-06-27 12:53:15 UTC (rev 6570)
@@ -20,6 +20,25 @@
import numpy as np
import numpy.fft
+# "large" composite numbers supported by FFTPACK
+LARGE_COMPOSITE_SIZES = [
+ 2**13,
+ 2**5 * 3**5,
+ 2**3 * 3**3 * 5**2,
+]
+SMALL_COMPOSITE_SIZES = [
+ 2,
+ 2*3*5,
+ 2*2*3*3,
+]
+# prime
+LARGE_PRIME_SIZES = [
+ 2011
+]
+SMALL_PRIME_SIZES = [
+ 29
+]
+
from numpy.random import rand
def random(size):
return rand(*size)
@@ -134,7 +153,6 @@
y = fftpack.zrfft(x)
assert_array_almost_equal(y,y2)
-
class TestDoubleFFT(_TestFFTBase):
def setUp(self):
self.cdt = np.cdouble
@@ -145,6 +163,10 @@
self.cdt = np.complex64
self.rdt = np.float32
+ @dec.knownfailureif(True, "single-precision FFT implementation is partially disabled, until accuracy issues with large prime powers are resolved")
+ def test_notice(self):
+ pass
+
class _TestIFFTBase(TestCase):
def test_definition(self):
x = np.array([1,2,3,4+1j,1,2,3,4+2j], self.cdt)
@@ -205,6 +227,31 @@
assert_array_almost_equal (y1, x)
assert_array_almost_equal (y2, x)
+ def test_size_accuracy(self):
+ # Sanity check for the accuracy for prime and non-prime sized inputs
+ if self.rdt == np.float32:
+ rtol = 1e-5
+ elif self.rdt == np.float64:
+ rtol = 1e-10
+
+ for size in LARGE_COMPOSITE_SIZES + LARGE_PRIME_SIZES:
+ np.random.seed(1234)
+ x = np.random.rand(size).astype(self.rdt)
+ y = ifft(fft(x))
+ self.failUnless(np.linalg.norm(x - y) < rtol*np.linalg.norm(x),
+ (size, self.rdt))
+ y = fft(ifft(x))
+ self.failUnless(np.linalg.norm(x - y) < rtol*np.linalg.norm(x),
+ (size, self.rdt))
+
+ x = (x + 1j*np.random.rand(size)).astype(self.cdt)
+ y = ifft(fft(x))
+ self.failUnless(np.linalg.norm(x - y) < rtol*np.linalg.norm(x),
+ (size, self.rdt))
+ y = fft(ifft(x))
+ self.failUnless(np.linalg.norm(x - y) < rtol*np.linalg.norm(x),
+ (size, self.rdt))
+
class TestDoubleIFFT(_TestIFFTBase):
def setUp(self):
self.cdt = np.cdouble
@@ -294,9 +341,28 @@
"Output dtype is %s, expected %s" % (y1.dtype, self.rdt))
self.failUnless(y2.dtype == self.rdt,
"Output dtype is %s, expected %s" % (y2.dtype, self.rdt))
- assert_array_almost_equal (y1, x, decimal=self.ndec)
- assert_array_almost_equal (y2, x, decimal=self.ndec)
+ assert_array_almost_equal (y1, x, decimal=self.ndec,
+ err_msg="size=%d" % size)
+ assert_array_almost_equal (y2, x, decimal=self.ndec,
+ err_msg="size=%d" % size)
+ def test_size_accuracy(self):
+ # Sanity check for the accuracy for prime and non-prime sized inputs
+ if self.rdt == np.float32:
+ rtol = 1e-5
+ elif self.rdt == np.float64:
+ rtol = 1e-10
+
+ for size in LARGE_COMPOSITE_SIZES + LARGE_PRIME_SIZES:
+ np.random.seed(1234)
+ x = np.random.rand(size).astype(self.rdt)
+ y = irfft(rfft(x))
+ self.failUnless(np.linalg.norm(x - y) < rtol*np.linalg.norm(x),
+ (size, self.rdt))
+ y = rfft(irfft(x))
+ self.failUnless(np.linalg.norm(x - y) < rtol*np.linalg.norm(x),
+ (size, self.rdt))
+
# self.ndec is bogus; we should have a assert_array_approx_equal for number of
# significant digits
class TestIRFFTDouble(_TestIRFFTBase):
@@ -331,6 +397,25 @@
y_r = np.array(fftn(x), np.complex64)
assert_array_almost_equal_nulp(y, y_r)
+ def test_size_accuracy(self):
+ for size in SMALL_COMPOSITE_SIZES + SMALL_PRIME_SIZES:
+ np.random.seed(1234)
+ x = np.random.rand(size, size) + 1j*np.random.rand(size, size)
+ y1 = fftn(x.astype(np.float32))
+ y2 = fftn(x.astype(np.float64)).astype(np.complex64)
+
+ self.failUnless(y1.dtype == np.complex64)
+ assert_array_almost_equal_nulp(y1, y2, 2000)
+
+ for size in LARGE_COMPOSITE_SIZES + LARGE_PRIME_SIZES:
+ np.random.seed(1234)
+ x = np.random.rand(size, 3) + 1j*np.random.rand(size, 3)
+ y1 = fftn(x.astype(np.float32))
+ y2 = fftn(x.astype(np.float64)).astype(np.complex64)
+
+ self.failUnless(y1.dtype == np.complex64)
+ assert_array_almost_equal_nulp(y1, y2, 2000)
+
class TestFftn(TestCase):
def test_definition(self):
@@ -513,7 +598,7 @@
class TestIfftnSingle(_TestIfftn):
dtype = np.float32
cdtype = np.complex64
- maxnlp = 2000
+ maxnlp = 3500
class TestLongDoubleFailure(TestCase):
def test_complex(self):
More information about the Scipy-svn
mailing list