[Scipy-svn] r6570 - in branches/0.8.x: doc/release scipy/fftpack scipy/fftpack/tests

scipy-svn@scip... scipy-svn@scip...
Sun Jun 27 07:53:15 CDT 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