[Scipy-svn] r6826 - in trunk/scipy/fftpack: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sun Sep 26 10:17:17 EDT 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