diff --git a/mkl_fft/__init__.py b/mkl_fft/__init__.py index 2d70e55..9b3f0a0 100644 --- a/mkl_fft/__init__.py +++ b/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'] diff --git a/mkl_fft/_pydfti.pyx b/mkl_fft/_pydfti.pyx index f594d4a..069170b 100644 --- a/mkl_fft/_pydfti.pyx +++ b/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.PyArray_EMPTY( - b_ndim, b_shape, b_type, 0) # 0 for C-contiguous + b_ndim, b_shape, 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.PyArray_EMPTY( - f_ndim, f_shape, f_type, 0) # 0 for C-contiguous + f_ndim, f_shape, 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 = 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_, 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 = 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_, 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_, 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_, axis_, f_arr) if (status): - raise ValueError("Internal error occurred") + raise ValueError("Internal error occurred, status={}".format(status)) return f_arr diff --git a/mkl_fft/setup.py b/mkl_fft/setup.py index c69fa6b..96ff036 100644 --- a/mkl_fft/setup.py +++ b/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', ] ) diff --git a/mkl_fft/tests/test_fftnd.py b/mkl_fft/tests/test_fftnd.py new file mode 100644 index 0000000..84f8781 --- /dev/null +++ b/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) diff --git a/setup.py b/setup.py index 3b191a3..f044b1e 100644 --- a/setup.py +++ b/setup.py @@ -28,7 +28,7 @@ MAJOR = 1 MINOR = 0 -MICRO = 2 +MICRO = 3 ISRELEASED = False VERSION = '%d.%d.%d' % (MAJOR, MINOR, MICRO)