Bug fixes

mkl_fft/__init__.py

@@ -27,5 +27,6 @@
 from ._pydfti import (fft, ifft, fft2, ifft2, fftn, ifftn, rfft, irfft,
                       rfft_numpy, irfft_numpy, rfftn_numpy, irfftn_numpy)
 
+__version__ = '1.0.3'
 __all__ = ['fft', 'ifft', 'fft2', 'ifft2', 'fftn', 'ifftn', 'rfft', 'irfft',
            'rfft_numpy', 'irfft_numpy', 'rfftn_numpy', 'irfftn_numpy']

mkl_fft/_pydfti.pyx

@@ -124,7 +124,7 @@ def ifft(x, n=None, axis=-1, overwrite_x=False):
 cdef cnp.ndarray pad_array(cnp.ndarray x_arr, cnp.npy_intp n, int axis, int realQ):
     "Internal utility to zero-pad input array along given axis"
     cdef cnp.ndarray b_arr "b_arrayObject"
-    cdef int x_type, b_type, b_ndim
+    cdef int x_type, b_type, b_ndim, x_arr_is_fortran
     cdef cnp.npy_intp *b_shape
 
     x_type = cnp.PyArray_TYPE(x_arr)
@@ -140,8 +140,9 @@ cdef cnp.ndarray pad_array(cnp.ndarray x_arr, cnp.npy_intp n, int axis, int real
     b_shape[axis] = n
 
     # allocating temporary buffer
+    x_arr_is_fortran = cnp.PyArray_CHKFLAGS(x_arr, cnp.NPY_F_CONTIGUOUS)
     b_arr = <cnp.ndarray> cnp.PyArray_EMPTY(
-        b_ndim, b_shape, <cnp.NPY_TYPES> b_type, 0) # 0 for C-contiguous
+        b_ndim, b_shape, <cnp.NPY_TYPES> b_type, x_arr_is_fortran) # 0 for C-contiguous
     PyMem_Free(b_shape)
 
     ind = [slice(0, None, None), ] * b_ndim
@@ -227,6 +228,7 @@ cdef cnp.ndarray __allocate_result(cnp.ndarray x_arr, long n_, long axis_, int f
     """
     cdef cnp.npy_intp *f_shape
     cdef cnp.ndarray f_arr "ff_arrayObject"
+    cdef int x_arr_is_fortran
 
     f_ndim = cnp.PyArray_NDIM(x_arr)
 
@@ -237,8 +239,9 @@ cdef cnp.ndarray __allocate_result(cnp.ndarray x_arr, long n_, long axis_, int f
         f_shape[axis_] = n_
 
     # allocating output buffer
+    x_arr_is_fortran = cnp.PyArray_CHKFLAGS(x_arr, cnp.NPY_F_CONTIGUOUS)
     f_arr = <cnp.ndarray> cnp.PyArray_EMPTY(
-        f_ndim, f_shape, <cnp.NPY_TYPES> f_type, 0) # 0 for C-contiguous
+        f_ndim, f_shape, <cnp.NPY_TYPES> f_type, x_arr_is_fortran) # 0 for C-contiguous
     PyMem_Free(f_shape);
 
     return f_arr
@@ -299,7 +302,7 @@ def _fft1d_impl(x, n=None, axis=-1, overwrite_arg=False, direction=+1):
             status = 1
 
         if status:
-            raise ValueError("Internal error")
+            raise ValueError("Internal error, status={}".format(status))
 
         n_max = <long> cnp.PyArray_DIM(x_arr, axis_)
         if (n_ < n_max):
@@ -346,7 +349,7 @@ def _fft1d_impl(x, n=None, axis=-1, overwrite_arg=False, direction=+1):
                         x_arr, n_, <int> axis_, f_arr)
 
         if (status):
-            raise ValueError("Internal error occurred")
+            raise ValueError("Internal error occurred, status={}".format(status))
 
         return f_arr
 
@@ -409,7 +412,7 @@ def _rrfft1d_impl(x, n=None, axis=-1, overwrite_arg=False, direction=+1):
             status = 1
 
         if status:
-            raise ValueError("Internal error")
+            raise ValueError("Internal error, status={}".format(status))
 
         n_max = <long> cnp.PyArray_DIM(x_arr, axis_)
         if (n_ < n_max):
@@ -434,7 +437,7 @@ def _rrfft1d_impl(x, n=None, axis=-1, overwrite_arg=False, direction=+1):
                status = float_float_mkl_rfft_out(x_arr, n_, <int> axis_, f_arr)
 
         if (status):
-            raise ValueError("Internal error occurred")
+            raise ValueError("Internal error occurred, status={}".format(status))
 
         return f_arr
 
@@ -488,7 +491,7 @@ def _rc_fft1d_impl(x, n=None, axis=-1, overwrite_arg=False):
         status = double_cdouble_mkl_fft1d_out(x_arr, n_, <int> axis_, f_arr, HALF_HARMONICS)
 
     if (status):
-        raise ValueError("Internal error occurred")
+        raise ValueError("Internal error occurred, with status={}".format(status))
 
     return f_arr
 
@@ -564,7 +567,7 @@ def _rc_ifft1d_impl(x, n=None, axis=-1, overwrite_arg=False):
             status = cdouble_double_mkl_irfft_out(x_arr, n_, <int> axis_, f_arr)
 
         if (status):
-            raise ValueError("Internal error occurred")
+            raise ValueError("Internal error occurred, status={}".format(status))
 
         return f_arr
 

mkl_fft/setup.py

@@ -31,7 +31,7 @@
 def configuration(parent_package='',top_path=None):
     from numpy.distutils.misc_util import Configuration
     from numpy.distutils.system_info import get_info
-    config = Configuration('mkl_fft', parent_package, top_path)
+    config = Configuration('', parent_package, top_path)
 
     pdir = dirname(__file__)
     wdir = join(pdir, 'src')
@@ -63,7 +63,7 @@ def configuration(parent_package='',top_path=None):
         libraries = libs,
         extra_compile_args = [
             '-DNDEBUG',
-            # '-g', '-O0', '-Wall', '-Wextra',
+            # '-g', '-O0', '-Wall', '-Wextra', '-DDEBUG',
         ]
     )
 

mkl_fft/tests/test_fftnd.py

@@ -0,0 +1,120 @@
+#!/usr/bin/env python
+# Copyright (c) 2017, Intel Corporation
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+#
+#     * Redistributions of source code must retain the above copyright notice,
+#       this list of conditions and the following disclaimer.
+#     * Redistributions in binary form must reproduce the above copyright
+#       notice, this list of conditions and the following disclaimer in the
+#       documentation and/or other materials provided with the distribution.
+#     * Neither the name of Intel Corporation nor the names of its contributors
+#       may be used to endorse or promote products derived from this software
+#       without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+from __future__ import division, absolute_import, print_function
+
+import numpy as np
+from numpy.testing import (
+        TestCase, run_module_suite, assert_, assert_raises, assert_equal,
+        assert_warns, assert_allclose)
+from numpy import random as rnd
+import sys
+import warnings
+
+import mkl_fft
+import numpy.fft.fftpack as np_fft
+
+reps_64 = (2**11)*np.finfo(np.float64).eps
+reps_32 = (2**11)*np.finfo(np.float32).eps
+atol_64 = (2**8)*np.finfo(np.float64).eps
+atol_32 = (2**8)*np.finfo(np.float32).eps
+
+def _get_rtol_atol(x):
+    dt = x.dtype
+    if dt == np.double or dt == np.complex128:
+        return reps_64, atol_64
+    elif dt == np.single or dt == np.complex64:
+        return reps_32, atol_32
+    else:
+        assert (dt == np.double or dt == np.complex128 or dt == np.single or dt == np.complex64), "Unexpected dtype {}".format(dt)
+        return reps_64, atol_64
+
+
+class Test_mklfft_matrix(TestCase):
+    def setUp(self):
+        rnd.seed(123456)
+        self.md = rnd.randn(256, 256)
+        self.mf = self.md.astype(np.float32)
+        self.mz = rnd.randn(256, 256*2).view(np.complex128)
+        self.mc = self.mz.astype(np.complex64)
+
+    def test_matrix1(self):
+        """fftn equals repeated fft"""
+        for ar in [self.md, self.mz, self.mf, self.mc]:
+            r_tol, a_tol = _get_rtol_atol(ar)
+            d = ar.copy()
+            t1 = mkl_fft.fftn(d)
+            t2 = mkl_fft.fft(mkl_fft.fft(d, axis=0), axis=1)
+            t3 = mkl_fft.fft(mkl_fft.fft(d, axis=1), axis=0)
+            assert_allclose(t1, t2, rtol=r_tol, atol=a_tol, err_msg = "failed test for dtype {}, max abs diff: {}".format(d.dtype, np.max(np.abs(t1-t2))))
+            assert_allclose(t1, t3, rtol=r_tol, atol=a_tol, err_msg = "failed test for dtype {}, max abs diff: {}".format(d.dtype, np.max(np.abs(t1-t3))))
+
+    def test_matrix2(self):
+        """ifftn(fftn(x)) is x"""
+        for ar in [self.md, self.mz, self.mf, self.mc]:
+            d = ar.copy()
+            r_tol, a_tol = _get_rtol_atol(d)
+            t = mkl_fft.ifftn(mkl_fft.fftn(d))
+            assert_allclose(d, t, rtol=r_tol, atol=a_tol, err_msg = "failed test for dtype {}, max abs diff: {}".format(d.dtype, np.max(np.abs(d-t))))
+
+    def test_matrix3(self):
+        """fftn(ifftn(x)) is x"""
+        for ar in [self.md, self.mz, self.mf, self.mc]:
+            d = ar.copy()
+            r_tol, a_tol = _get_rtol_atol(d)
+            t = mkl_fft.fftn(mkl_fft.ifftn(d))
+            assert_allclose(d, t, rtol=r_tol, atol=a_tol, err_msg = "failed test for dtype {}, max abs diff: {}".format(d.dtype, np.max(np.abs(d-t))))
+
+
+    def test_matrix4(self):
+        """fftn of strided array is same as fftn of a contiguous copy"""
+        for ar in [self.md, self.mz, self.mf, self.mc]:
+            r_tol, a_tol = _get_rtol_atol(ar)
+            d_strided = ar[::2,::2]
+            d_contig = d_strided.copy()
+            t_strided = mkl_fft.fftn(d_strided)
+            t_contig = mkl_fft.fftn(d_contig)
+            assert_allclose(t_strided, t_contig, rtol=r_tol, atol=a_tol)
+
+
+class Test_Regressions(TestCase):
+
+    def setUp(self):
+        rnd.seed(123456)
+        self.ad = rnd.randn(32, 17, 23)
+        self.af = self.ad.astype(np.float32)
+        self.az = rnd.randn(32, 17, 23*2).view(np.complex128)
+        self.ac = self.az.astype(np.complex64)
+
+    def test_cf_contig(self):
+        """fft of F-contiguous array is the same as of C-contiguous with same data"""
+        for ar in [self.ad, self.af, self.az, self.ac]:
+            r_tol, a_tol = _get_rtol_atol(ar)
+            d_ccont = ar.copy()
+            d_fcont = np.asfortranarray(d_ccont)
+            f1 = mkl_fft.fft(d_ccont)
+            f2 = mkl_fft.fft(d_fcont)
+            assert_allclose(f1, f2, rtol=r_tol, atol=a_tol)

setup.py

@@ -28,7 +28,7 @@
 
 MAJOR = 1
 MINOR = 0
-MICRO = 2
+MICRO = 3
 ISRELEASED = False
 
 VERSION = '%d.%d.%d' % (MAJOR, MINOR, MICRO)