[Scipy-svn] r6826 - in trunk/scipy/fftpack: . tests

scipy-svn@scip... scipy-svn@scip...
Sun Sep 26 09:17:17 CDT 2010


Author: ptvirtan
Date: 2010-09-26 09:17:17 -0500 (Sun, 26 Sep 2010)
New Revision: 6826

Modified:
   trunk/scipy/fftpack/basic.py
   trunk/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.

(cherry picked from r6570)

Modified: trunk/scipy/fftpack/basic.py
===================================================================
--- trunk/scipy/fftpack/basic.py	2010-09-26 02:51:25 UTC (rev 6825)
+++ trunk/scipy/fftpack/basic.py	2010-09-26 14:17:17 UTC (rev 6826)
@@ -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: trunk/scipy/fftpack/tests/test_basic.py
===================================================================
--- trunk/scipy/fftpack/tests/test_basic.py	2010-09-26 02:51:25 UTC (rev 6825)
+++ trunk/scipy/fftpack/tests/test_basic.py	2010-09-26 14:17:17 UTC (rev 6826)
@@ -12,7 +12,8 @@
 """
 
 from numpy.testing import assert_, assert_equal, assert_array_almost_equal, \
-        assert_array_almost_equal_nulp, assert_raises, run_module_suite, TestCase
+        assert_array_almost_equal_nulp, assert_raises, run_module_suite, \
+        TestCase, dec
 from scipy.fftpack import ifft,fft,fftn,ifftn,rfft,irfft, fft2
 from scipy.fftpack import _fftpack as fftpack
 
@@ -21,6 +22,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)
@@ -136,7 +156,6 @@
             y = fftpack.zrfft(x)
             assert_array_almost_equal(y,y2)
 
-
 class TestDoubleFFT(_TestFFTBase):
     def setUp(self):
         self.cdt = np.cdouble
@@ -147,6 +166,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 setUp(self):
         np.random.seed(1234)
@@ -210,6 +233,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
@@ -303,9 +351,28 @@
                     "Output dtype is %s, expected %s" % (y1.dtype, self.rdt))
             self.assertTrue(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):
@@ -346,6 +413,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 setUp(self):
         np.random.seed(1234)
@@ -529,7 +615,7 @@
 class TestIfftnSingle(_TestIfftn):
     dtype = np.float32
     cdtype = np.complex64
-    maxnlp = 2000
+    maxnlp = 3500
 
 class TestLongDoubleFailure(TestCase):
     def setUp(self):



More information about the Scipy-svn mailing list