Skip to content

GitLab

Commit 395a0c89 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

change normalization parameters

parent d77cdeb5
Hide whitespace changes
Inline Side-by-side
Showing with 168 additions and 125 deletions
bench_nd.py
View file @ 395a0c89
......@@ -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))
......
pocketfft_hdronly.h
View file @ 395a0c89
......@@ -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 <cmath>
......@@ -48,9 +48,6 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#include <memory>
#include <vector>
#include <complex>
#if defined(_WIN32)
#include <malloc.h>
#endif
#ifdef POCKETFFT_OPENMP
#include <omp.h>
#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<typename T> 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<T *>(res);
void *ptr = malloc(num*sizeof(T)+64);
if (!ptr) throw bad_alloc();
T *res = reinterpret_cast<T *>
((reinterpret_cast<size_t>(ptr) & ~(size_t(63))) + 64);
(reinterpret_cast<void**>(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<T *>(res);
}
static void dealloc(T *ptr)
{ free(ptr); }
{ if (ptr) free((reinterpret_cast<void**>(ptr))[-1]); }
#endif
public:
......
pypocketfft.cc
View file @ 395a0c89
......@@ -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<double> args; \
if (dtype.is(T2)) return func<float> args; \
if (dtype.is(T3)) return func<long double> args; \
if (dtype.is(T3)) return func<ldbl_t> args; \
throw runtime_error("unsupported data type");
template<typename T> 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<typename T> 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<T>(inorm, N);
}
template<typename T> 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<complex<T>>(dims);
......@@ -86,28 +108,29 @@ template<typename T> py::array xfftn_internal(const py::array &in,
auto d_out=reinterpret_cast<complex<T> *>(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<T>(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<typename T> 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<typename T> py::array rfftn_internal(const py::array &in,
auto d_out=reinterpret_cast<complex<T> *>(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<T>(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<typename T> 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<T>(dims);
......@@ -141,26 +165,27 @@ template<typename T> py::array xrfft_scipy(const py::array &in,
auto d_out=reinterpret_cast<T *>(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<T>(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<typename T> 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<typename T> py::array irfftn_internal(const py::array &in,
auto d_out=reinterpret_cast<T *>(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<T>(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<typename T> 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<T>(dims);
......@@ -200,15 +226,16 @@ template<typename T> py::array hartley_internal(const py::array &in,
auto d_out=reinterpret_cast<T *>(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<T>(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 @@ template<typename T>py::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);
......
stress.py
View file @ 395a0c89
......@@ -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
......
test.py
View file @ 395a0c89
......@@ -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)