Commit c1cddc18 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

improve Hartley documentation and names

parent 15e37d0f
...@@ -16,10 +16,6 @@ ...@@ -16,10 +16,6 @@
#include "pocketfft_hdronly.h" #include "pocketfft_hdronly.h"
//
// Python interface
//
namespace { namespace {
using namespace std; using namespace std;
...@@ -175,7 +171,8 @@ template<typename T> py::array r2c_internal(const py::array &in, ...@@ -175,7 +171,8 @@ template<typename T> py::array r2c_internal(const py::array &in,
{ {
py::gil_scoped_release release; py::gil_scoped_release release;
T fct = norm_fct<T>(inorm, dims_in, axes); T fct = norm_fct<T>(inorm, dims_in, axes);
pocketfft::r2c(dims_in, s_in, s_out, axes, forward, d_in, d_out, fct, nthreads); pocketfft::r2c(dims_in, s_in, s_out, axes, forward, d_in, d_out, fct,
nthreads);
} }
return res; return res;
} }
...@@ -201,8 +198,8 @@ template<typename T> py::array r2r_fftpack_internal(const py::array &in, ...@@ -201,8 +198,8 @@ template<typename T> py::array r2r_fftpack_internal(const py::array &in,
{ {
py::gil_scoped_release release; py::gil_scoped_release release;
T fct = norm_fct<T>(inorm, dims, axes); T fct = norm_fct<T>(inorm, dims, axes);
pocketfft::r2r_fftpack(dims, s_in, s_out, axes, !input_halfcomplex, forward, d_in, d_out, pocketfft::r2r_fftpack(dims, s_in, s_out, axes, !input_halfcomplex, forward,
fct, nthreads); d_in, d_out, fct, nthreads);
} }
return res; return res;
} }
...@@ -233,7 +230,8 @@ template<typename T> py::array c2r_internal(const py::array &in, ...@@ -233,7 +230,8 @@ template<typename T> py::array c2r_internal(const py::array &in,
{ {
py::gil_scoped_release release; py::gil_scoped_release release;
T fct = norm_fct<T>(inorm, dims_out, axes); T fct = norm_fct<T>(inorm, dims_out, axes);
pocketfft::c2r(dims_out, s_in, s_out, axes, forward, d_in, d_out, fct, nthreads); pocketfft::c2r(dims_out, s_in, s_out, axes, forward, d_in, d_out, fct,
nthreads);
} }
return res; return res;
} }
...@@ -245,7 +243,7 @@ py::array c2r(const py::array &in, py::object axes_, size_t lastsize, ...@@ -245,7 +243,7 @@ py::array c2r(const py::array &in, py::object axes_, size_t lastsize,
inorm, nthreads)) inorm, nthreads))
} }
template<typename T> py::array hartley_internal(const py::array &in, template<typename T> py::array separable_hartley_internal(const py::array &in,
py::object axes_, int inorm, bool inplace, size_t nthreads) py::object axes_, int inorm, bool inplace, size_t nthreads)
{ {
auto dims(copy_shape(in)); auto dims(copy_shape(in));
...@@ -263,11 +261,11 @@ template<typename T> py::array hartley_internal(const py::array &in, ...@@ -263,11 +261,11 @@ template<typename T> py::array hartley_internal(const py::array &in,
return res; return res;
} }
py::array hartley(const py::array &in, py::object axes_, int inorm, py::array separable_hartley(const py::array &in, py::object axes_, int inorm,
bool inplace, size_t nthreads) bool inplace, size_t nthreads)
{ {
DISPATCH(in, f64, f32, flong, hartley_internal, (in, axes_, inorm, inplace, DISPATCH(in, f64, f32, flong, separable_hartley_internal, (in, axes_, inorm,
nthreads)) inplace, nthreads))
} }
template<typename T>py::array complex2hartley(const py::array &in, template<typename T>py::array complex2hartley(const py::array &in,
...@@ -295,17 +293,11 @@ template<typename T>py::array complex2hartley(const py::array &in, ...@@ -295,17 +293,11 @@ template<typename T>py::array complex2hartley(const py::array &in,
return out; return out;
} }
py::array mycomplex2hartley(const py::array &in, py::array genuine_hartley(const py::array &in, py::object axes_, int inorm,
const py::array &tmp, py::object axes_, bool inplace)
{
DISPATCH(in, f64, f32, flong, complex2hartley, (in, tmp, axes_, inplace))
}
py::array hartley2(const py::array &in, py::object axes_, int inorm,
bool inplace, size_t nthreads) bool inplace, size_t nthreads)
{ {
return mycomplex2hartley(in, r2c(in, axes_, true, inorm, nthreads), axes_, auto tmp = r2c(in, axes_, true, inorm, nthreads);
inplace); DISPATCH(in, f64, f32, flong, complex2hartley, (in, tmp, axes_, inplace))
} }
const char *pypocketfft_DS = R"""(Fast Fourier and Hartley transforms. const char *pypocketfft_DS = R"""(Fast Fourier and Hartley transforms.
...@@ -446,10 +438,42 @@ np.ndarray (same shape and data type as `a`) ...@@ -446,10 +438,42 @@ np.ndarray (same shape and data type as `a`)
The transformed data. The shape is identical to that of the input array. The transformed data. The shape is identical to that of the input array.
)"""; )""";
const char *hartley_DS = R"""(Performs a Hartley transform. const char *separable_hartley_DS = R"""(Performs a separable Hartley transform.
For every requested axis, a 1D forward Fourier transform is carried out, For every requested axis, a 1D forward Fourier transform is carried out, and
and the sum of real and imaginary parts of the result is stored in the output the real and imaginary parts of the result are added before the next axis is
array. processed.
Parameters
----------
a : numpy.ndarray (np.float32 or np.float64)
The input data
axes : list of integers
The axes along which the transform is carried out.
If not set, all axes will be transformed.
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.
nthreads : int
Number of threads to use. If 0, use the system default (typically governed
by the `OMP_NUM_THREADS` environment variable).
Returns
-------
np.ndarray (same shape and data type as `a`)
The transformed data
)""";
const char *genuine_hartley_DS = R"""(Performs a full Hartley transform.
A full Fourier transform is carried out over the requested axes, and the
sum of real and imaginary parts of the result is stored in the output
array. For a single transformed axis, this is identical to `separable_hartley`,
but when transforming multiple axes, the results are different.
Parameters Parameters
---------- ----------
...@@ -493,10 +517,8 @@ PYBIND11_MODULE(pypocketfft, m) ...@@ -493,10 +517,8 @@ PYBIND11_MODULE(pypocketfft, m)
m.def("r2r_fftpack",&r2r_fftpack, r2r_fftpack_DS, "a"_a, "axes"_a, m.def("r2r_fftpack",&r2r_fftpack, r2r_fftpack_DS, "a"_a, "axes"_a,
"input_halfcomplex"_a, "forward"_a, "inorm"_a=0, "inplace"_a=false, "input_halfcomplex"_a, "forward"_a, "inorm"_a=0, "inplace"_a=false,
"nthreads"_a=1); "nthreads"_a=1);
m.def("hartley",&hartley, hartley_DS, "a"_a, "axes"_a=py::none(), "inorm"_a=0, m.def("separable_hartley",&separable_hartley, separable_hartley_DS, "a"_a,
"inplace"_a=false, "nthreads"_a=1); "axes"_a=py::none(), "inorm"_a=0, "inplace"_a=false, "nthreads"_a=1);
m.def("hartley2",&hartley2, "a"_a, "axes"_a=py::none(), "inorm"_a=0, m.def("genuine_hartley",&genuine_hartley, genuine_hartley_DS, "a"_a,
"inplace"_a=false, "nthreads"_a=1); "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);
} }
...@@ -77,22 +77,22 @@ def test(): ...@@ -77,22 +77,22 @@ def test():
if err > fmaxerrf: if err > fmaxerrf:
fmaxerrf = err fmaxerrf = err
print("fmaxerrf:", fmaxerrf, shape, axes) print("fmaxerrf:", fmaxerrf, shape, axes)
b=pypocketfft.hartley(pypocketfft.hartley(a.real,axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads) b=pypocketfft.separable_hartley(pypocketfft.separable_hartley(a.real,axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads)
err = _l2error(a.real,b) err = _l2error(a.real,b)
if err > hmaxerr: if err > hmaxerr:
hmaxerr = err hmaxerr = err
print("hmaxerr:", hmaxerr, shape, axes) print("hmaxerr:", hmaxerr, shape, axes)
b=pypocketfft.hartley2(pypocketfft.hartley2(a.real,axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads) b=pypocketfft.genuine_hartley(pypocketfft.genuine_hartley(a.real,axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads)
err = _l2error(a.real,b) err = _l2error(a.real,b)
if err > hmaxerr: if err > hmaxerr:
hmaxerr = err hmaxerr = err
print("hmaxerr:", hmaxerr, shape, axes) print("hmaxerr:", hmaxerr, shape, axes)
b=pypocketfft.hartley(pypocketfft.hartley(a.real.astype(np.float32),axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads) b=pypocketfft.separable_hartley(pypocketfft.separable_hartley(a.real.astype(np.float32),axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads)
err = _l2error(a.real.astype(np.float32),b) err = _l2error(a.real.astype(np.float32),b)
if err > hmaxerrf: if err > hmaxerrf:
hmaxerrf = err hmaxerrf = err
print("hmaxerrf:", hmaxerrf, shape, axes) print("hmaxerrf:", hmaxerrf, shape, axes)
b=pypocketfft.hartley2(pypocketfft.hartley2(a.real.astype(np.float32),axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads) b=pypocketfft.genuine_hartley(pypocketfft.genuine_hartley(a.real.astype(np.float32),axes=axes,nthreads=nthreads),axes=axes,inorm=2,nthreads=nthreads)
err = _l2error(a.real.astype(np.float32),b) err = _l2error(a.real.astype(np.float32),b)
if err > hmaxerrf: if err > hmaxerrf:
hmaxerrf = err hmaxerrf = err
......
...@@ -127,9 +127,9 @@ def test_identity_r2(shp): ...@@ -127,9 +127,9 @@ def test_identity_r2(shp):
assert_(_l2error(rfftn(irfftn(a),inorm=2), a)<1e-15) assert_(_l2error(rfftn(irfftn(a),inorm=2), a)<1e-15)
@pmp("shp", shapes2D+shapes3D) @pmp("shp", shapes2D+shapes3D)
def test_hartley2(shp): def test_genuine_hartley(shp):
a=np.random.rand(*shp)-0.5 a=np.random.rand(*shp)-0.5
v1 = pypocketfft.hartley2(a) v1 = pypocketfft.genuine_hartley(a)
v2 = fftn(a.astype(np.complex128)) v2 = fftn(a.astype(np.complex128))
v2 = v2.real+v2.imag v2 = v2.real+v2.imag
assert_(_l2error(v1, v2)<1e-15) assert_(_l2error(v1, v2)<1e-15)
...@@ -137,20 +137,20 @@ def test_hartley2(shp): ...@@ -137,20 +137,20 @@ def test_hartley2(shp):
@pmp("shp", shapes) @pmp("shp", shapes)
def test_hartley_identity(shp): def test_hartley_identity(shp):
a=np.random.rand(*shp)-0.5 a=np.random.rand(*shp)-0.5
v1 = pypocketfft.hartley(pypocketfft.hartley(a))/a.size v1 = pypocketfft.separable_hartley(pypocketfft.separable_hartley(a))/a.size
assert_(_l2error(a, v1)<1e-15) assert_(_l2error(a, v1)<1e-15)
@pmp("shp", shapes) @pmp("shp", shapes)
def test_hartley2_identity(shp): def test_genuine_hartley_identity(shp):
a=np.random.rand(*shp)-0.5 a=np.random.rand(*shp)-0.5
v1 = pypocketfft.hartley2(pypocketfft.hartley2(a))/a.size v1 = pypocketfft.genuine_hartley(pypocketfft.genuine_hartley(a))/a.size
assert_(_l2error(a, v1)<1e-15) assert_(_l2error(a, v1)<1e-15)
v1 = a.copy() v1 = a.copy()
assert_ (pypocketfft.hartley2(pypocketfft.hartley2(v1, inplace=True), inorm=2, inplace=True) is v1) assert_ (pypocketfft.genuine_hartley(pypocketfft.genuine_hartley(v1, inplace=True), inorm=2, inplace=True) is v1)
assert_(_l2error(a, v1)<1e-15) assert_(_l2error(a, v1)<1e-15)
@pmp("shp", shapes2D) @pmp("shp", shapes2D)
@pmp("axes", ((0,),(1,),(0,1),(1,0))) @pmp("axes", ((0,),(1,),(0,1),(1,0)))
def test_hartley2_2D(shp, axes): def test_genuine_hartley_2D(shp, axes):
a=np.random.rand(*shp)-0.5 a=np.random.rand(*shp)-0.5
assert_(_l2error(pypocketfft.hartley2(pypocketfft.hartley2(a, axes=axes), axes=axes, inorm=2),a)<1e-15) assert_(_l2error(pypocketfft.genuine_hartley(pypocketfft.genuine_hartley(a, axes=axes), axes=axes, inorm=2),a)<1e-15)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment