Added FFT tests to CI and make test.

ci/scripts/test.sh

@@ -52,6 +52,7 @@ python -c 'import hysop; print hysop'
 
 python "$HYSOP_DIR/core/arrays/tests/test_array.py"
 $HYSOP_DIR/fields/tests/test_cartesian.sh
+python "$HYSOP_DIR/numerics/tests/test_fft.py"
 python "$HYSOP_DIR/operator/tests/test_analytic.py"
 python "$HYSOP_DIR/operator/tests/test_transpose.py"
 python "$HYSOP_DIR/operator/tests/test_derivative.py"

cmake/hysop_tests.cmake

@@ -52,7 +52,7 @@ endmacro()
 set(testDir ${HYSOP_BUILD_PYTHONPATH})
 
 # === Set the list of all directories which may contain tests ===
-set(py_src_dirs core/arrays fields operator)
+set(py_src_dirs core/arrays numerics fields operator)
 
 # === Create the files list from all directories in py_src_dirs ===
 

hysop/numerics/fft/gpyfft_fft.py

@@ -4,7 +4,7 @@ FFT iterface for fast Fourier Transforms using CLFFT backend.
 :class:`~hysop.numerics.GpyFFTPlan`
 """
 
-import warnings, struct
+import warnings, struct, primefac
 import numpy as np
 from abc import abstractmethod
 
@@ -250,6 +250,7 @@ class GpyFFTPlan(FFTPlanI):
                     (out_array.strides[t_axes_in[0]] == out_array.dtype.itemsize)), \
                     'Inplace real transforms need stride 1 for first transform axis.'
         
+        self.check_transform_shape(t_shape)
         plan = GFFT.create_plan(self.context, t_shape[::-1])
         plan.inplace       = t_inplace
         plan.strides_in    = t_strides_in[::-1]
@@ -342,7 +343,7 @@ t_inplace, t_shape, t_precision,
 layout_in, layout_out, plan.scale_forward, plan.scale_backward, 
 self.pre_callback_src, self.post_callback_src)
             print msg
-    
+
     def set_callbacks(self, plan, axes, N,
             in_array, out_array, fake_input, fake_output,
             layout_in, layout_out, **kwds):
@@ -380,6 +381,24 @@ self.pre_callback_src, self.post_callback_src)
         if (post_src is not None):
             plan.set_callback('post_callback', post_src, 'post', user_data=user_data)
         return (in_data, out_data)
+
+    @classmethod
+    def check_transform_shape(self, shape):
+        """Check that clFFT can handle the logical transform size."""
+        valid_factors = {2,3,5,7,11,13}
+        for Ni in shape:
+            factors = tuple( primefac.primefac(int(Ni)) )
+            invalid_factors = set(factors) - valid_factors
+            if invalid_factors:
+                factorization = ' * '.join('{}^{}'.format(factor, factors.count(factor)) 
+                                                                for factor in set(factors))
+                candidates = ', '.join(str(vf) for vf in valid_factors)
+                msg ='\nInvalid transform shape {} for clFFT:'
+                msg+='\n   {} = {}'
+                msg+='\nOnly {} prime factors are available.'
+                msg+='\n'
+                msg=msg.format(shape, Ni, factorization, candidates)
+                raise ValueError(msg)
     
     @classmethod
     def calculate_transform_strides(cls, taxes, array):

hysop/numerics/tests/test_fft.py

@@ -123,7 +123,7 @@ class TestFFT(object):
         print '\n FORWARD R2C: real to hermitian complex transform'
         for (shape, cshape, rshape, N, Nc, Nr) in self.iter_shapes():
             print '   RFFT shape={:9s} '.format(str(shape)+':'),
-            Href = np.random.rand(N).astype(dtype).reshape(shape)
+            Href = np.random.rand(*shape).astype(dtype).reshape(shape)
             results = {}
             for (kind, implementations) in self.implementations.iteritems():
                 for (name, impl) in implementations.iteritems():
@@ -444,7 +444,7 @@ class TestFFT(object):
         print '\n R2C-C2R transform'
         for (shape, cshape, rshape, N, Nc, Nr) in self.iter_shapes():
             print '   X - IRFFT(RFFT(X)) shape={:9s} '.format(str(shape)+':'),
-            Href = np.random.rand(N).astype(dtype).reshape(shape)
+            Href = np.random.rand(*shape).astype(dtype).reshape(shape)
             results = {}
             for (kind, implementations) in self.implementations.iteritems():
                 for (name, impl) in implementations.iteritems():
@@ -467,8 +467,8 @@ class TestFFT(object):
                     results[name] = np.max(np.abs(Href-H0))
             check_distances(results)
         
-        return
         print '\n R2R-R2R transforms'
+
         types = ['I  ','II ','III','IV ']
         for (itype,stype) in enumerate(types, 1):
             print '\n DCT-{}: real to real discrete cosine transform {}'.format(
@@ -476,64 +476,81 @@ class TestFFT(object):
             ttype = 'COS{}'.format(itype)
             for (shape, cshape, rshape, N, Nc, Nr) in self.iter_shapes():
                 print '   X - I{}({}(X)) shape={:9s} '.format(ttype, ttype, str(shape)+':'),
-                Href = np.random.rand(N).astype(dtype).reshape(shape)
-                H0   = np.empty(shape=shape,  dtype=dtype)
-                H1   = np.empty(shape=shape,  dtype=dtype)
-                H0, H1, Href = self.simd_align(H0, H1, Href)
+                if (itype==1): # real size is 2*(N-1)
+                    shape = mk_shape(shape, -1, shape[-1] + 1)
+                Href = np.random.rand(*shape).astype(dtype).reshape(shape)
                 results = {}
-                for (name, impl) in self.implementations.iteritems():
-                    iitype = [1,3,2,4][itype-1]
-                    if dtype not in impl.supported_ftypes:
-                        continue
-                    if itype not in impl.supported_cosine_transforms:
-                        continue
-                    if iitype not in impl.supported_cosine_transforms:
-                        continue
-                    forward  = impl.dct(a=H0, out=H1, type=itype)
-                    backward = impl.idct(a=H1, out=H0, type=itype)
-                    H0[...] = Href
-                    forward.execute()
-                    backward.execute()
-                    results[name] = np.max(np.abs(Href-H0))
+                for (kind, implementations) in self.implementations.iteritems():
+                    for (name, impl) in implementations.iteritems():
+                        iitype = [1,3,2,4][itype-1]
+                        if dtype not in impl.supported_ftypes:
+                            continue
+                        if itype not in impl.supported_cosine_transforms:
+                            continue
+                        if iitype not in impl.supported_cosine_transforms:
+                            continue
+                        D0 = impl.backend.empty(shape=shape, dtype=dtype)
+                        D1 = impl.backend.empty(shape=shape, dtype=dtype)
+                        forward  = impl.dct(a=D0, out=D1, type=itype).setup()
+                        backward = impl.idct(a=D1, out=D0, type=itype).setup()
+                        assert forward.input_array   is D0
+                        assert forward.output_array  is D1
+                        assert backward.input_array  is D1
+                        assert backward.output_array is D0
+                        D0[...] = Href
+                        D1[...] = np.nan
+                        forward.execute()
+                        D0[...] = np.nan
+                        backward.execute()
+                        H0 = D0.get()
+                        results[name] = np.max(np.abs(Href-H0))
                 check_distances(results)
         
         for (itype,stype) in enumerate(types, 1):
-            print '\n DST-{}: real to real discrete sine transform {}'.format(
+            print '\n DST-{}: real to real discrete sinine transform {}'.format(
                     stype.strip(), itype)
             ttype = 'SIN{}'.format(itype)
             for (shape, cshape, rshape, N, Nc, Nr) in self.iter_shapes():
                 print '   X - I{}({}(X)) shape={:9s} '.format(ttype, ttype, str(shape)+':'),
-                Href = np.random.rand(N).astype(dtype).reshape(shape)
-                H0   = np.empty(shape=shape,  dtype=dtype)
-                H1   = np.empty(shape=shape,  dtype=dtype)
-                H0, H1, Href = self.simd_align(H0, H1, Href)
+                if (itype==1): # real size is 2*(N+1)
+                    shape = mk_shape(shape, -1, shape[-1] - 1)
+                Href = np.random.rand(*shape).astype(dtype).reshape(shape)
                 results = {}
-                for (name, impl) in self.implementations.iteritems():
-                    iitype = [1,3,2,4][itype-1]
-                    if dtype not in impl.supported_ftypes:
-                        continue
-                    if itype not in impl.supported_sine_transforms:
-                        continue
-                    if iitype not in impl.supported_sine_transforms:
-                        continue
-                    forward  = impl.dst(a=H0, out=H1, type=itype)
-                    backward = impl.idst(a=H1, out=H0, type=itype)
-                    H0[...] = Href
-                    forward.execute()
-                    backward.execute()
-                    results[name] = np.max(np.abs(Href-H0))
+                for (kind, implementations) in self.implementations.iteritems():
+                    for (name, impl) in implementations.iteritems():
+                        iitype = [1,3,2,4][itype-1]
+                        if dtype not in impl.supported_ftypes:
+                            continue
+                        if itype not in impl.supported_sine_transforms:
+                            continue
+                        if iitype not in impl.supported_sine_transforms:
+                            continue
+                        D0 = impl.backend.empty(shape=shape, dtype=dtype)
+                        D1 = impl.backend.empty(shape=shape, dtype=dtype)
+                        forward  = impl.dst(a=D0, out=D1, type=itype).setup()
+                        backward = impl.idst(a=D1, out=D0, type=itype).setup()
+                        assert forward.input_array   is D0
+                        assert forward.output_array  is D1
+                        assert backward.input_array  is D1
+                        assert backward.output_array is D0
+                        D0[...] = Href
+                        D1[...] = np.nan
+                        forward.execute()
+                        D0[...] = np.nan
+                        backward.execute()
+                        H0 = D0.get()
+                        results[name] = np.max(np.abs(Href-H0))
                 check_distances(results)
     
     def iter_shapes(self):
         minj=(2,2)
-        maxj=(3,3)
-        # maxj=(6,6)
+        maxj=(5,5)
         msg = ('EVEN', 'ODD')
         for i in xrange(2):
             base = 2+i
             print '  '+msg[i]
             for j1 in xrange(minj[i],maxj[i]):
-                shape = (2,4,3,base**j1,)
+                shape = (3,2,base**j1,)
                 cshape = list(shape)
                 cshape[-1] = cshape[-1]//2 + 1
                 cshape = tuple(cshape)
@@ -566,15 +583,6 @@ class TestFFT(object):
             ss += (s,)
             failed = (Eeps>report_eps)
             if failed:
-                print
-                print 'r0={}, r1={}'.format(r0, r1)
-                print
-                print results[r0]
-                print
-                print results[r1]
-                print
-                import sys
-                sys.exit(1)
                 shape=results[r0].shape
                 failures.setdefault(tag, []).append((r0, r1, shape, Einf, Eeps))
                 
@@ -585,12 +593,6 @@ class TestFFT(object):
             msg='Some implementations did not agree on the result !'
             raise RuntimeError(msg)
 
-    def simd_align(self, *arrays):
-        aligned_arrays = ()
-        for array in arrays:
-            aligned_arrays += (pyfftw.byte_align(array),)
-        return aligned_arrays
-
     def report_failures(self, failures):
         print
         print '== TEST FAILURES REPORT =='
@@ -635,9 +637,9 @@ class TestFFT(object):
         # not testing np.longdouble because only one implementation supports it
         # ie. we cannot compare results between different implementations
         failures = {}
-        for dtype in (np.float64,):
+        for dtype in (np.float32, np.float64,):
             self._test_forward_backward_1d(dtype=dtype)
-            # self._test_1d(dtype=dtype, failures=failures.setdefault(dtype.__name__, {}))
+            self._test_1d(dtype=dtype, failures=failures.setdefault(dtype.__name__, {}))
         self.report_failures(failures)
     
 if __name__ == '__main__':