diff --git a/bench_nd.py b/bench_nd.py index f532b0e207307034dbee59e81dee6a2f0fb03b9b..ab391751ddd601eb5872a5cee4cf68e4450e1575 100644 --- a/bench_nd.py +++ b/bench_nd.py @@ -5,7 +5,7 @@ import pypocketfft from time import time import matplotlib.pyplot as plt -nthreads=0 +nthreads=1 def _l2error(a,b): return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2)) @@ -26,7 +26,7 @@ def bench_nd_fftn(ndim, nmax, ntry, tp, nrepeat, filename=""): b=pypocketfft.fftn(a,nthreads=nthreads) t1=time() tmin_pp = min(tmin_pp,t1-t0) - a2=pypocketfft.ifftn(b,fct=1./a.size) + a2=pypocketfft.ifftn(b,inorm=2) assert(_l2error(a,a2)<(2.5e-15 if tp=='c16' else 6e-7)) res.append(tmin_pp/tmin_np) plt.title("t(pypocketfft / numpy 1.17), {}D, {}, max_extent={}".format(ndim, str(tp), nmax)) diff --git a/pocketfft_hdronly.h b/pocketfft_hdronly.h index b467629c15307e15c3cd459a6f170d01af33bd8d..372ad2d41a05245476ab868be676f50fd9cff9cc 100644 --- a/pocketfft_hdronly.h +++ b/pocketfft_hdronly.h @@ -34,11 +34,11 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #define POCKETFFT_HDRONLY_H #ifndef __cplusplus -#error This file is C++ and requires a C++ compiler +#error This file is C++ and requires a C++ compiler. #endif #if !(__cplusplus >= 201103L || _MSVC_LANG+0L >= 201103L) -#error This file requires at least C++11 support +#error This file requires at least C++11 support. #endif #include @@ -48,9 +48,6 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include #include #include -#if defined(_WIN32) -#include -#endif #ifdef POCKETFFT_OPENMP #include #endif @@ -83,7 +80,7 @@ constexpr bool FORWARD = true, #ifndef POCKETFFT_NO_VECTORS #define POCKETFFT_NO_VECTORS #if defined(__INTEL_COMPILER) -// do nothing. This is necessary because this compiler also sets __GNUC__ +// do nothing. This is necessary because this compiler also sets __GNUC__. #elif defined(__clang__) #if __clang__>=5 #undef POCKETFFT_NO_VECTORS @@ -141,27 +138,19 @@ template class arr } static void dealloc(T *ptr) { free(ptr); } -#elif defined(_WIN32) +#else // portable emulation static T *ralloc(size_t num) { if (num==0) return nullptr; - void *res = _aligned_malloc(num*sizeof(T), 64); - if (!res) throw bad_alloc(); - return reinterpret_cast(res); + void *ptr = malloc(num*sizeof(T)+64); + if (!ptr) throw bad_alloc(); + T *res = reinterpret_cast + ((reinterpret_cast(ptr) & ~(size_t(63))) + 64); + (reinterpret_cast(res))[-1] = ptr; + return res; } static void dealloc(T *ptr) - { _aligned_free(ptr); } -#else - static T *ralloc(size_t num) - { - if (num==0) return nullptr; - void *res(nullptr); - if (posix_memalign(&res, 64, num*sizeof(T))!=0) - throw bad_alloc(); - return reinterpret_cast(res); - } - static void dealloc(T *ptr) - { free(ptr); } + { if (ptr) free((reinterpret_cast(ptr))[-1]); } #endif public: diff --git a/pypocketfft.cc b/pypocketfft.cc index 6fb321563b61ec91efef90b958f521ad430cbb88..fabe4a27d8e71652af70c674afe3b2284cfde83f 100644 --- a/pypocketfft.cc +++ b/pypocketfft.cc @@ -27,12 +27,16 @@ using namespace pocketfft; namespace py = pybind11; +// Only instantiate long double transforms if they offer more precision +using ldbl_t = typename std::conditional< + sizeof(long double)==sizeof(double), double, long double>::type; + auto c64 = py::dtype("complex64"); auto c128 = py::dtype("complex128"); -auto c256 = py::dtype("complex256"); +auto clong = py::dtype("longcomplex"); auto f32 = py::dtype("float32"); auto f64 = py::dtype("float64"); -auto f128 = py::dtype("float128"); +auto flong = py::dtype("longfloat"); shape_t copy_shape(const py::array &arr) { @@ -72,11 +76,29 @@ shape_t makeaxes(const py::array &in, py::object axes) auto dtype = arr.dtype(); \ if (dtype.is(T1)) return func args; \ if (dtype.is(T2)) return func args; \ - if (dtype.is(T3)) return func args; \ + if (dtype.is(T3)) return func args; \ throw runtime_error("unsupported data type"); +template T norm_fct(int inorm, size_t N) + { + if (inorm==0) return T(1); + if (inorm==2) return T(1/ldbl_t(N)); + if (inorm==1) return T(1/sqrt(ldbl_t(N))); + throw invalid_argument("invalid value for inorm (must be 0, 1, or 2)"); + } + +template T norm_fct(int inorm, const shape_t &shape, + const shape_t &axes) + { + if (inorm==0) return T(1); + size_t N(1); + for (auto a: axes) + N *= shape[a]; + return norm_fct(inorm, N); + } + template py::array xfftn_internal(const py::array &in, - const shape_t &axes, long double fct, bool inplace, bool fwd, size_t nthreads) + const shape_t &axes, int inorm, bool inplace, bool fwd, size_t nthreads) { auto dims(copy_shape(in)); py::array res = inplace ? in : py::array_t>(dims); @@ -86,28 +108,29 @@ template py::array xfftn_internal(const py::array &in, auto d_out=reinterpret_cast *>(res.mutable_data()); { py::gil_scoped_release release; - c2c(dims, s_in, s_out, axes, fwd, d_in, d_out, T(fct), nthreads); + T fct = norm_fct(inorm, dims, axes); + c2c(dims, s_in, s_out, axes, fwd, d_in, d_out, fct, nthreads); } return res; } -py::array xfftn(const py::array &a, py::object axes, double fct, bool inplace, - bool fwd, size_t nthreads) +py::array xfftn(const py::array &a, py::object axes, int inorm, + bool inplace, bool fwd, size_t nthreads) { - DISPATCH(a, c128, c64, c256, xfftn_internal, (a, makeaxes(a, axes), fct, + DISPATCH(a, c128, c64, clong, xfftn_internal, (a, makeaxes(a, axes), inorm, inplace, fwd, nthreads)) } -py::array fftn(const py::array &a, py::object axes, double fct, bool inplace, - size_t nthreads) - { return xfftn(a, axes, fct, inplace, true, nthreads); } +py::array fftn(const py::array &a, py::object axes, int inorm, + bool inplace, size_t nthreads) + { return xfftn(a, axes, inorm, inplace, true, nthreads); } -py::array ifftn(const py::array &a, py::object axes, double fct, bool inplace, - size_t nthreads) - { return xfftn(a, axes, fct, inplace, false, nthreads); } +py::array ifftn(const py::array &a, py::object axes, int inorm, + bool inplace, size_t nthreads) + { return xfftn(a, axes, inorm, inplace, false, nthreads); } template py::array rfftn_internal(const py::array &in, - py::object axes_, long double fct, size_t nthreads) + py::object axes_, int inorm, size_t nthreads) { auto axes = makeaxes(in, axes_); auto dims_in(copy_shape(in)), dims_out(dims_in); @@ -119,19 +142,20 @@ template py::array rfftn_internal(const py::array &in, auto d_out=reinterpret_cast *>(res.mutable_data()); { py::gil_scoped_release release; - r2c(dims_in, s_in, s_out, axes, d_in, d_out, T(fct), nthreads); + T fct = norm_fct(inorm, dims_in, axes); + r2c(dims_in, s_in, s_out, axes, d_in, d_out, fct, nthreads); } return res; } -py::array rfftn(const py::array &in, py::object axes_, double fct, +py::array rfftn(const py::array &in, py::object axes_, int inorm, size_t nthreads) { - DISPATCH(in, f64, f32, f128, rfftn_internal, (in, axes_, fct, nthreads)) + DISPATCH(in, f64, f32, flong, rfftn_internal, (in, axes_, inorm, nthreads)) } template py::array xrfft_scipy(const py::array &in, - size_t axis, long double fct, bool inplace, bool fwd, size_t nthreads) + size_t axis, int inorm, bool inplace, bool fwd, size_t nthreads) { auto dims(copy_shape(in)); py::array res = inplace ? in : py::array_t(dims); @@ -141,26 +165,27 @@ template py::array xrfft_scipy(const py::array &in, auto d_out=reinterpret_cast(res.mutable_data()); { py::gil_scoped_release release; - r2r_fftpack(dims, s_in, s_out, axis, fwd, d_in, d_out, T(fct), nthreads); + T fct = norm_fct(inorm, dims[axis]); + r2r_fftpack(dims, s_in, s_out, axis, fwd, d_in, d_out, fct, nthreads); } return res; } -py::array rfft_scipy(const py::array &in, size_t axis, double fct, bool inplace, - size_t nthreads) +py::array rfft_scipy(const py::array &in, size_t axis, int inorm, + bool inplace, size_t nthreads) { - DISPATCH(in, f64, f32, f128, xrfft_scipy, (in, axis, fct, inplace, true, + DISPATCH(in, f64, f32, flong, xrfft_scipy, (in, axis, inorm, inplace, true, nthreads)) } -py::array irfft_scipy(const py::array &in, size_t axis, double fct, +py::array irfft_scipy(const py::array &in, size_t axis, int inorm, bool inplace, size_t nthreads) { - DISPATCH(in, f64, f32, f128, xrfft_scipy, (in, axis, fct, inplace, false, + DISPATCH(in, f64, f32, flong, xrfft_scipy, (in, axis, inorm, inplace, false, nthreads)) } template py::array irfftn_internal(const py::array &in, - py::object axes_, size_t lastsize, long double fct, size_t nthreads) + py::object axes_, size_t lastsize, int inorm, size_t nthreads) { auto axes = makeaxes(in, axes_); size_t axis = axes.back(); @@ -176,20 +201,21 @@ template py::array irfftn_internal(const py::array &in, auto d_out=reinterpret_cast(res.mutable_data()); { py::gil_scoped_release release; - c2r(dims_out, s_in, s_out, axes, d_in, d_out, T(fct), nthreads); + T fct = norm_fct(inorm, dims_out, axes); + c2r(dims_out, s_in, s_out, axes, d_in, d_out, fct, nthreads); } return res; } py::array irfftn(const py::array &in, py::object axes_, size_t lastsize, - double fct, size_t nthreads) + int inorm, size_t nthreads) { - DISPATCH(in, c128, c64, c256, irfftn_internal, (in, axes_, lastsize, fct, + DISPATCH(in, c128, c64, clong, irfftn_internal, (in, axes_, lastsize, inorm, nthreads)) } template py::array hartley_internal(const py::array &in, - py::object axes_, long double fct, bool inplace, size_t nthreads) + py::object axes_, int inorm, bool inplace, size_t nthreads) { auto dims(copy_shape(in)); py::array res = inplace ? in : py::array_t(dims); @@ -200,15 +226,16 @@ template py::array hartley_internal(const py::array &in, auto d_out=reinterpret_cast(res.mutable_data()); { py::gil_scoped_release release; - r2r_hartley(dims, s_in, s_out, axes, d_in, d_out, T(fct), nthreads); + T fct = norm_fct(inorm, dims, axes); + r2r_hartley(dims, s_in, s_out, axes, d_in, d_out, fct, nthreads); } return res; } -py::array hartley(const py::array &in, py::object axes_, double fct, +py::array hartley(const py::array &in, py::object axes_, int inorm, bool inplace, size_t nthreads) { - DISPATCH(in, f64, f32, f128, hartley_internal, (in, axes_, fct, inplace, + DISPATCH(in, f64, f32, flong, hartley_internal, (in, axes_, inorm, inplace, nthreads)) } @@ -261,13 +288,13 @@ templatepy::array complex2hartley(const py::array &in, py::array mycomplex2hartley(const py::array &in, const py::array &tmp, py::object axes_, bool inplace) { - DISPATCH(in, f64, f32, f128, complex2hartley, (in, tmp, axes_, inplace)) + DISPATCH(in, f64, f32, flong, complex2hartley, (in, tmp, axes_, inplace)) } -py::array hartley2(const py::array &in, py::object axes_, double fct, +py::array hartley2(const py::array &in, py::object axes_, int inorm, bool inplace, size_t nthreads) { - return mycomplex2hartley(in, rfftn(in, axes_, fct, nthreads), axes_, + return mycomplex2hartley(in, rfftn(in, axes_, inorm, nthreads), axes_, inplace); } @@ -293,8 +320,12 @@ a : numpy.ndarray (np.complex64 or np.complex128) axes : list of integers The axes along which the FFT is carried out. If not set, all axes will be transformed. -fct : float - Normalization factor +inorm : int + Normalization type + 0 : no normalization + 1 : divide by sqrt(N) + 2 : divide by N + where N is the product of the lengths of the transformed axes. inplace : bool if False, returns the result in a new array and leaves the input unchanged. if True, stores the result in the input array and returns a handle to it. @@ -317,8 +348,12 @@ a : numpy.ndarray (np.complex64 or np.complex128) axes : list of integers The axes along which the FFT is carried out. If not set, all axes will be transformed. -fct : float - Normalization factor +inorm : int + Normalization type + 0 : no normalization + 1 : divide by sqrt(N) + 2 : divide by N + where N is the product of the lengths of the transformed axes. inplace : bool if False, returns the result in a new array and leaves the input unchanged. if True, stores the result in the input array and returns a handle to it. @@ -341,8 +376,12 @@ a : numpy.ndarray (np.float32 or np.float64) axes : list of integers The axes along which the FFT is carried out. If not set, all axes will be transformed in ascending order. -fct : float - Normalization factor +inorm : int + Normalization type + 0 : no normalization + 1 : divide by sqrt(N) + 2 : divide by N + where N is the product of the lengths of the transformed input axes. nthreads : int Number of threads to use. If 0, use the system default (typically governed by the `OMP_NUM_THREADS` environment variable). @@ -363,8 +402,12 @@ a : numpy.ndarray (np.float32 or np.float64) The input data axis : int The axis along which the FFT is carried out. -fct : float - Normalization factor +inorm : int + Normalization type + 0 : no normalization + 1 : divide by sqrt(N) + 2 : divide by N + where N is the length of `axis`. inplace : bool if False, returns the result in a new array and leaves the input unchanged. if True, stores the result in the input array and returns a handle to it. @@ -391,8 +434,12 @@ axes : list of integers If not set, all axes will be transformed in ascending order. lastsize : the output size of the last axis to be transformed. If the corresponding input axis has size n, this can be 2*n-2 or 2*n-1. -fct : float - Normalization factor +inorm : int + Normalization type + 0 : no normalization + 1 : divide by sqrt(N) + 2 : divide by N + where N is the product of the lengths of the transformed output axes. nthreads : int Number of threads to use. If 0, use the system default (typically governed by the `OMP_NUM_THREADS` environment variable). @@ -414,8 +461,12 @@ a : numpy.ndarray (np.float32 or np.float64) FFTPACK half-complex order, i.e. `a[0].re, a[1].re, a[1].im, a[2].re ...`. axis : int The axis along which the FFT is carried out. -fct : float - Normalization factor +inorm : int + Normalization type + 0 : no normalization + 1 : divide by sqrt(N) + 2 : divide by N + where N is the length of `axis`. inplace : bool if False, returns the result in a new array and leaves the input unchanged. if True, stores the result in the input array and returns a handle to it. @@ -441,8 +492,12 @@ a : numpy.ndarray (np.float32 or np.float64) axes : list of integers The axes along which the transform is carried out. If not set, all axes will be transformed. -fct : float - Normalization factor +inorm : int + Normalization type + 0 : no normalization + 1 : divide by sqrt(N) + 2 : divide by N + where N is the product of the lengths of the transformed axes. inplace : bool if False, returns the result in a new array and leaves the input unchanged. if True, stores the result in the input array and returns a handle to it. @@ -463,21 +518,21 @@ PYBIND11_MODULE(pypocketfft, m) using namespace pybind11::literals; m.doc() = pypocketfft_DS; - m.def("fftn",&fftn, fftn_DS, "a"_a, "axes"_a=py::none(), "fct"_a=1., + m.def("fftn",&fftn, fftn_DS, "a"_a, "axes"_a=py::none(), "inorm"_a=0, "inplace"_a=false, "nthreads"_a=1); - m.def("ifftn",&ifftn, ifftn_DS, "a"_a, "axes"_a=py::none(), "fct"_a=1., + m.def("ifftn",&ifftn, ifftn_DS, "a"_a, "axes"_a=py::none(), "inorm"_a=0, "inplace"_a=false, "nthreads"_a=1); - m.def("rfftn",&rfftn, rfftn_DS, "a"_a, "axes"_a=py::none(), "fct"_a=1., + m.def("rfftn",&rfftn, rfftn_DS, "a"_a, "axes"_a=py::none(), "inorm"_a=0, "nthreads"_a=1); - m.def("rfft_scipy",&rfft_scipy, rfft_scipy_DS, "a"_a, "axis"_a, "fct"_a=1., + m.def("rfft_scipy",&rfft_scipy, rfft_scipy_DS, "a"_a, "axis"_a, "inorm"_a=0, "inplace"_a=false, "nthreads"_a=1); m.def("irfftn",&irfftn, irfftn_DS, "a"_a, "axes"_a=py::none(), "lastsize"_a=0, - "fct"_a=1., "nthreads"_a=1); - m.def("irfft_scipy",&irfft_scipy, irfft_scipy_DS, "a"_a, "axis"_a, "fct"_a=1., - "inplace"_a=false, "nthreads"_a=1); - m.def("hartley",&hartley, hartley_DS, "a"_a, "axes"_a=py::none(), "fct"_a=1., + "inorm"_a=0, "nthreads"_a=1); + m.def("irfft_scipy",&irfft_scipy, irfft_scipy_DS, "a"_a, "axis"_a, + "inorm"_a=0, "inplace"_a=false, "nthreads"_a=1); + m.def("hartley",&hartley, hartley_DS, "a"_a, "axes"_a=py::none(), "inorm"_a=0, "inplace"_a=false, "nthreads"_a=1); - m.def("hartley2",&hartley2, "a"_a, "axes"_a=py::none(), "fct"_a=1., + m.def("hartley2",&hartley2, "a"_a, "axes"_a=py::none(), "inorm"_a=0, "inplace"_a=false, "nthreads"_a=1); m.def("complex2hartley",&mycomplex2hartley, "in"_a, "tmp"_a, "axes"_a, "inplace"_a=false); diff --git a/stress.py b/stress.py index 422c4a45ba1265af42e9e335ae07b563049c11cc..66644c31aa36a39f4e9973b9cd5199ca83976955 100644 --- a/stress.py +++ b/stress.py @@ -21,34 +21,33 @@ def test(): nax = np.random.randint(1,ndim+1) axes = axes[:nax] lastsize = shape[axes[-1]] - fct = 1./np.prod(np.take(shape, axes)) a=np.random.rand(*shape)-0.5 + 1j*np.random.rand(*shape)-0.5j - b=pypocketfft.ifftn(pypocketfft.fftn(a,axes=axes,nthreads=nthreads),axes=axes,fct=fct,nthreads=nthreads) + b=pypocketfft.ifftn(pypocketfft.fftn(a,axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads) err = _l2error(a,b) if err > cmaxerr: cmaxerr = err print("cmaxerr:", cmaxerr, shape, axes) - b=pypocketfft.irfftn(pypocketfft.rfftn(a.real,axes=axes,nthreads=nthreads),axes=axes,fct=fct,lastsize=lastsize,nthreads=nthreads) + b=pypocketfft.irfftn(pypocketfft.rfftn(a.real,axes=axes,nthreads=nthreads),axes=axes,inorm=2,lastsize=lastsize,nthreads=nthreads) err = _l2error(a.real,b) if err > fmaxerr: fmaxerr = err print("fmaxerr:", fmaxerr, shape, axes) - b=pypocketfft.ifftn(pypocketfft.fftn(a.astype(np.complex64),axes=axes,nthreads=nthreads),axes=axes,fct=fct,nthreads=nthreads) + b=pypocketfft.ifftn(pypocketfft.fftn(a.astype(np.complex64),axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads) err = _l2error(a.astype(np.complex64),b) if err > cmaxerrf: cmaxerrf = err print("cmaxerrf:", cmaxerrf, shape, axes) - b=pypocketfft.irfftn(pypocketfft.rfftn(a.real.astype(np.float32),axes=axes,nthreads=nthreads),axes=axes,fct=fct,lastsize=lastsize,nthreads=nthreads) + b=pypocketfft.irfftn(pypocketfft.rfftn(a.real.astype(np.float32),axes=axes,nthreads=nthreads),axes=axes,inorm=2,lastsize=lastsize,nthreads=nthreads) err = _l2error(a.real.astype(np.float32),b) if err > fmaxerrf: fmaxerrf = err print("fmaxerrf:", fmaxerrf, shape, axes) - b=pypocketfft.hartley(pypocketfft.hartley(a.real,axes=axes,nthreads=nthreads),axes=axes,fct=fct,nthreads=nthreads) + b=pypocketfft.hartley(pypocketfft.hartley(a.real,axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads) err = _l2error(a.real,b) if err > hmaxerr: hmaxerr = err print("hmaxerr:", hmaxerr, shape, axes) - b=pypocketfft.hartley(pypocketfft.hartley(a.real.astype(np.float32),axes=axes,nthreads=nthreads),axes=axes,fct=fct,nthreads=nthreads) + b=pypocketfft.hartley(pypocketfft.hartley(a.real.astype(np.float32),axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads) err = _l2error(a.real.astype(np.float32),b) if err > hmaxerrf: hmaxerrf = err diff --git a/test.py b/test.py index 080756d80aa6e1115e4c1e693b296ed743c8f559..d0f547cabea8737e98b6554243a0d902a74fee36 100644 --- a/test.py +++ b/test.py @@ -2,6 +2,7 @@ import pypocketfft import pyfftw import numpy as np import pytest +from numpy.testing import assert_ pmp = pytest.mark.parametrize @@ -15,69 +16,70 @@ def _l2error(a,b): return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2)) @pmp("len", len1D) -def test1D(len): +@pmp("inorm", [0,1,2]) +def test1D(len, inorm): a=np.random.rand(len)-0.5 + 1j*np.random.rand(len)-0.5j b=a.astype(np.complex64) c=a.astype(np.complex256) - assert(_l2error(a, pypocketfft.ifftn(pypocketfft.fftn(c))*np.float128(1.)/len)<1e-18) - assert(_l2error(a, pypocketfft.ifftn(pypocketfft.fftn(a),fct=1./len))<1.5e-15) - assert(_l2error(a.real, pypocketfft.irfftn(pypocketfft.rfftn(a.real),fct=1./len,lastsize=len))<1.5e-15) + assert_(_l2error(a, pypocketfft.ifftn(pypocketfft.fftn(c,inorm=inorm), inorm=2-inorm))<1e-18) + assert_(_l2error(a, pypocketfft.ifftn(pypocketfft.fftn(a,inorm=inorm), inorm=2-inorm))<1.5e-15) + assert_(_l2error(a.real, pypocketfft.irfftn(pypocketfft.rfftn(a.real,inorm=inorm), inorm=2-inorm,lastsize=len))<1.5e-15) tmp=a.copy() - assert (pypocketfft.ifftn(pypocketfft.fftn(tmp, inplace=True), fct=1./len, inplace=True) is tmp) - assert(_l2error(tmp, a)<1.5e-15) - assert(_l2error(b, pypocketfft.ifftn(pypocketfft.fftn(b),fct=1./len))<6e-7) - assert(_l2error(b.real, pypocketfft.irfftn(pypocketfft.rfftn(b.real),fct=1./len,lastsize=len))<6e-7) + assert_ (pypocketfft.ifftn(pypocketfft.fftn(tmp, inplace=True, inorm=inorm), inplace=True, inorm=2-inorm) is tmp) + assert_(_l2error(tmp, a)<1.5e-15) + assert_(_l2error(b, pypocketfft.ifftn(pypocketfft.fftn(b, inorm=inorm), inorm=2-inorm))<6e-7) + assert_(_l2error(b.real, pypocketfft.irfftn(pypocketfft.rfftn(b.real, inorm=inorm), lastsize=len, inorm=2-inorm))<6e-7) tmp=b.copy() - assert (pypocketfft.ifftn(pypocketfft.fftn(tmp, inplace=True), fct=1./len, inplace=True) is tmp) - assert(_l2error(tmp, b)<6e-7) + assert_ (pypocketfft.ifftn(pypocketfft.fftn(tmp, inplace=True, inorm=inorm), inplace=True, inorm=2-inorm) is tmp) + assert_(_l2error(tmp, b)<6e-7) @pmp("shp", shapes) @pmp("nthreads", (0,1,2)) def test_fftn(shp, nthreads): a=np.random.rand(*shp)-0.5 + 1j*np.random.rand(*shp)-0.5j - assert(_l2error(pyfftw.interfaces.numpy_fft.fftn(a), pypocketfft.fftn(a, nthreads=nthreads))<1e-15) + assert_(_l2error(pyfftw.interfaces.numpy_fft.fftn(a), pypocketfft.fftn(a, nthreads=nthreads))<1e-15) a=a.astype(np.complex64) - assert(_l2error(pyfftw.interfaces.numpy_fft.fftn(a), pypocketfft.fftn(a, nthreads=nthreads))<5e-7) + assert_(_l2error(pyfftw.interfaces.numpy_fft.fftn(a), pypocketfft.fftn(a, nthreads=nthreads))<5e-7) @pmp("shp", shapes2D) @pmp("axes", ((0,),(1,),(0,1),(1,0))) def test_fftn2D(shp, axes): a=np.random.rand(*shp)-0.5 + 1j*np.random.rand(*shp)-0.5j - assert(_l2error(pyfftw.interfaces.numpy_fft.fftn(a, axes=axes), pypocketfft.fftn(a, axes=axes))<1e-15) + assert_(_l2error(pyfftw.interfaces.numpy_fft.fftn(a, axes=axes), pypocketfft.fftn(a, axes=axes))<1e-15) a=a.astype(np.complex64) - assert(_l2error(pyfftw.interfaces.numpy_fft.fftn(a, axes=axes), pypocketfft.fftn(a, axes=axes))<5e-7) + assert_(_l2error(pyfftw.interfaces.numpy_fft.fftn(a, axes=axes), pypocketfft.fftn(a, axes=axes))<5e-7) @pmp("shp", shapes) def test_rfftn(shp): a=np.random.rand(*shp)-0.5 - assert(_l2error(pyfftw.interfaces.numpy_fft.rfftn(a), pypocketfft.rfftn(a))<1e-15) + assert_(_l2error(pyfftw.interfaces.numpy_fft.rfftn(a), pypocketfft.rfftn(a))<1e-15) a=a.astype(np.float32) - assert(_l2error(pyfftw.interfaces.numpy_fft.rfftn(a), pypocketfft.rfftn(a))<5e-7) + assert_(_l2error(pyfftw.interfaces.numpy_fft.rfftn(a), pypocketfft.rfftn(a))<5e-7) @pmp("shp", shapes) def test_rfft_scipy(shp): for i in range(len(shp)): a=np.random.rand(*shp)-0.5 - assert(_l2error(pyfftw.interfaces.scipy_fftpack.rfft(a,axis=i), pypocketfft.rfft_scipy(a,axis=i))<1e-15) - assert(_l2error(pyfftw.interfaces.scipy_fftpack.irfft(a,axis=i), pypocketfft.irfft_scipy(a,axis=i,fct=1./a.shape[i]))<1e-15) + assert_(_l2error(pyfftw.interfaces.scipy_fftpack.rfft(a,axis=i), pypocketfft.rfft_scipy(a,axis=i))<1e-15) + assert_(_l2error(pyfftw.interfaces.scipy_fftpack.irfft(a,axis=i), pypocketfft.irfft_scipy(a,axis=i,inorm=2))<1e-15) @pmp("shp", shapes2D) @pmp("axes", ((0,),(1,),(0,1),(1,0))) def test_rfftn2D(shp, axes): a=np.random.rand(*shp)-0.5 - assert(_l2error(pyfftw.interfaces.numpy_fft.rfftn(a, axes=axes), pypocketfft.rfftn(a, axes=axes))<1e-15) + assert_(_l2error(pyfftw.interfaces.numpy_fft.rfftn(a, axes=axes), pypocketfft.rfftn(a, axes=axes))<1e-15) a=a.astype(np.float32) - assert(_l2error(pyfftw.interfaces.numpy_fft.rfftn(a, axes=axes), pypocketfft.rfftn(a, axes=axes))<5e-7) + assert_(_l2error(pyfftw.interfaces.numpy_fft.rfftn(a, axes=axes), pypocketfft.rfftn(a, axes=axes))<5e-7) @pmp("shp", shapes) def test_identity(shp): a=np.random.rand(*shp)-0.5 + 1j*np.random.rand(*shp)-0.5j - assert(_l2error(pypocketfft.ifftn(pypocketfft.fftn(a),fct=1./a.size), a)<1.5e-15) + assert_(_l2error(pypocketfft.ifftn(pypocketfft.fftn(a),inorm=2), a)<1.5e-15) tmp=a.copy() - assert (pypocketfft.ifftn(pypocketfft.fftn(tmp, inplace=True), fct=1./a.size, inplace=True) is tmp) - assert(_l2error(tmp, a)<1.5e-15) + assert_ (pypocketfft.ifftn(pypocketfft.fftn(tmp, inplace=True), inorm=2, inplace=True) is tmp) + assert_(_l2error(tmp, a)<1.5e-15) a=a.astype(np.complex64) - assert(_l2error(pypocketfft.ifftn(pypocketfft.fftn(a),fct=1./a.size), a)<6e-7) + assert_(_l2error(pypocketfft.ifftn(pypocketfft.fftn(a),inorm=2), a)<6e-7) @pmp("shp", shapes) def test_identity_r(shp): @@ -85,15 +87,14 @@ def test_identity_r(shp): b=a.astype(np.float32) for ax in range(a.ndim): n = a.shape[ax] - assert(_l2error(pypocketfft.irfftn(pypocketfft.rfftn(a,(ax,)),(ax,),lastsize=n,fct=1./n), a)<1e-15) - assert(_l2error(pypocketfft.irfftn(pypocketfft.rfftn(b,(ax,)),(ax,),lastsize=n,fct=1./n), b)<5e-7) + assert_(_l2error(pypocketfft.irfftn(pypocketfft.rfftn(a,(ax,)),(ax,),lastsize=n,inorm=2), a)<1e-15) + assert_(_l2error(pypocketfft.irfftn(pypocketfft.rfftn(b,(ax,)),(ax,),lastsize=n,inorm=2), b)<5e-7) @pmp("shp", shapes) def test_identity_r2(shp): a=np.random.rand(*shp)-0.5 + 1j*np.random.rand(*shp)-0.5j a=pypocketfft.rfftn(pypocketfft.irfftn(a)) - fct = 1./pypocketfft.irfftn(a).size - assert(_l2error(pypocketfft.rfftn(pypocketfft.irfftn(a),fct=fct), a)<1e-15) + assert_(_l2error(pypocketfft.rfftn(pypocketfft.irfftn(a),inorm=2), a)<1e-15) @pmp("shp", shapes2D+shapes3D) def test_hartley2(shp): @@ -101,26 +102,25 @@ def test_hartley2(shp): v1 = pypocketfft.hartley2(a) v2 = pypocketfft.fftn(a.astype(np.complex128)) v2 = v2.real+v2.imag - assert(_l2error(v1, v2)<1e-15) + assert_(_l2error(v1, v2)<1e-15) @pmp("shp", shapes) def test_hartley_identity(shp): a=np.random.rand(*shp)-0.5 v1 = pypocketfft.hartley(pypocketfft.hartley(a))/a.size - assert(_l2error(a, v1)<1e-15) + assert_(_l2error(a, v1)<1e-15) @pmp("shp", shapes) def test_hartley2_identity(shp): a=np.random.rand(*shp)-0.5 v1 = pypocketfft.hartley2(pypocketfft.hartley2(a))/a.size - assert(_l2error(a, v1)<1e-15) + assert_(_l2error(a, v1)<1e-15) v1 = a.copy() - assert (pypocketfft.hartley2(pypocketfft.hartley2(v1, inplace=True), fct=1./a.size, inplace=True) is v1) - assert(_l2error(a, v1)<1e-15) + assert_ (pypocketfft.hartley2(pypocketfft.hartley2(v1, inplace=True), inorm=2, inplace=True) is v1) + assert_(_l2error(a, v1)<1e-15) @pmp("shp", shapes2D) @pmp("axes", ((0,),(1,),(0,1),(1,0))) def test_hartley2_2D(shp, axes): a=np.random.rand(*shp)-0.5 - fct = 1./np.prod(np.take(shp,axes)) - assert(_l2error(pypocketfft.hartley2(pypocketfft.hartley2(a, axes=axes), axes=axes, fct=fct),a)<1e-15) + assert_(_l2error(pypocketfft.hartley2(pypocketfft.hartley2(a, axes=axes), axes=axes, inorm=2),a)<1e-15)