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

towards reproducible test runs

parent 69b0aed4
...@@ -71,7 +71,8 @@ dtypes = [np.float32, np.float64, ...@@ -71,7 +71,8 @@ dtypes = [np.float32, np.float64,
@pmp("inorm", [0, 1, 2]) @pmp("inorm", [0, 1, 2])
@pmp("dtype", dtypes) @pmp("dtype", dtypes)
def test1D(len, inorm, dtype): def test1D(len, inorm, dtype):
a = np.random.rand(len)-0.5 + 1j*np.random.rand(len)-0.5j rng = np.random.default_rng(np.random.SeedSequence(42))
a = rng.random(len)-0.5 + 1j*rng.random(len)-0.5j
a = a.astype(ctype[dtype]) a = a.astype(ctype[dtype])
eps = tol[dtype] eps = tol[dtype]
assert_(_l2error(a, ifftn(fftn(a, inorm=inorm), inorm=2-inorm)) < eps) assert_(_l2error(a, ifftn(fftn(a, inorm=inorm), inorm=2-inorm)) < eps)
...@@ -91,7 +92,8 @@ def test1D(len, inorm, dtype): ...@@ -91,7 +92,8 @@ def test1D(len, inorm, dtype):
@pmp("nthreads", (0, 1, 2)) @pmp("nthreads", (0, 1, 2))
@pmp("inorm", [0, 1, 2]) @pmp("inorm", [0, 1, 2])
def test_fftn(shp, nthreads, inorm): def test_fftn(shp, nthreads, inorm):
a = np.random.rand(*shp)-0.5 + 1j*np.random.rand(*shp)-0.5j rng = np.random.default_rng(np.random.SeedSequence(42))
a = rng.random(shp)-0.5 + 1j*rng.random(shp)-0.5j
assert_(_l2error(a, ifftn(fftn(a, nthreads=nthreads, inorm=inorm), assert_(_l2error(a, ifftn(fftn(a, nthreads=nthreads, inorm=inorm),
nthreads=nthreads, inorm=2-inorm)) < 1e-15) nthreads=nthreads, inorm=2-inorm)) < 1e-15)
a = a.astype(np.complex64) a = a.astype(np.complex64)
...@@ -103,7 +105,8 @@ def test_fftn(shp, nthreads, inorm): ...@@ -103,7 +105,8 @@ def test_fftn(shp, nthreads, inorm):
@pmp("axes", ((0,), (1,), (0, 1), (1, 0))) @pmp("axes", ((0,), (1,), (0, 1), (1, 0)))
@pmp("inorm", [0, 1, 2]) @pmp("inorm", [0, 1, 2])
def test_fftn2D(shp, axes, inorm): def test_fftn2D(shp, axes, inorm):
a = np.random.rand(*shp)-0.5 + 1j*np.random.rand(*shp)-0.5j rng = np.random.default_rng(np.random.SeedSequence(42))
a = rng.random(shp)-0.5 + 1j*rng.random(shp)-0.5j
assert_(_l2error(a, ifftn(fftn(a, axes=axes, inorm=inorm), assert_(_l2error(a, ifftn(fftn(a, axes=axes, inorm=inorm),
axes=axes, inorm=2-inorm)) < 1e-15) axes=axes, inorm=2-inorm)) < 1e-15)
a = a.astype(np.complex64) a = a.astype(np.complex64)
...@@ -113,7 +116,8 @@ def test_fftn2D(shp, axes, inorm): ...@@ -113,7 +116,8 @@ def test_fftn2D(shp, axes, inorm):
@pmp("shp", shapes) @pmp("shp", shapes)
def test_rfftn(shp): def test_rfftn(shp):
a = np.random.rand(*shp)-0.5 rng = np.random.default_rng(np.random.SeedSequence(42))
a = rng.random(shp)-0.5
tmp1 = rfftn(a) tmp1 = rfftn(a)
tmp2 = fftn(a) tmp2 = fftn(a)
part = tuple(slice(0,tmp1.shape[i]) for i in range(tmp1.ndim)) part = tuple(slice(0,tmp1.shape[i]) for i in range(tmp1.ndim))
...@@ -128,7 +132,7 @@ def test_rfftn(shp): ...@@ -128,7 +132,7 @@ def test_rfftn(shp):
# @pmp("shp", shapes) # @pmp("shp", shapes)
# def test_rfft_scipy(shp): # def test_rfft_scipy(shp):
# for i in range(len(shp)): # for i in range(len(shp)):
# a = np.random.rand(*shp)-0.5 # a = rng.random(shp)-0.5
# assert_(_l2error(pyfftw.interfaces.scipy_fftpack.rfft(a, axis=i), # assert_(_l2error(pyfftw.interfaces.scipy_fftpack.rfft(a, axis=i),
# rfft_scipy(a, axis=i)) < 1e-15) # rfft_scipy(a, axis=i)) < 1e-15)
# assert_(_l2error(pyfftw.interfaces.scipy_fftpack.irfft(a, axis=i), # assert_(_l2error(pyfftw.interfaces.scipy_fftpack.irfft(a, axis=i),
...@@ -138,7 +142,8 @@ def test_rfftn(shp): ...@@ -138,7 +142,8 @@ def test_rfftn(shp):
@pmp("shp", shapes2D) @pmp("shp", shapes2D)
@pmp("axes", ((0,), (1,), (0, 1), (1, 0))) @pmp("axes", ((0,), (1,), (0, 1), (1, 0)))
def test_rfftn2D(shp, axes): def test_rfftn2D(shp, axes):
a = np.random.rand(*shp)-0.5 rng = np.random.default_rng(np.random.SeedSequence(42))
a = rng.random(shp)-0.5
tmp1 = rfftn(a,axes=axes) tmp1 = rfftn(a,axes=axes)
tmp2 = fftn(a,axes=axes) tmp2 = fftn(a,axes=axes)
part = tuple(slice(0,tmp1.shape[i]) for i in range(tmp1.ndim)) part = tuple(slice(0,tmp1.shape[i]) for i in range(tmp1.ndim))
...@@ -152,7 +157,8 @@ def test_rfftn2D(shp, axes): ...@@ -152,7 +157,8 @@ def test_rfftn2D(shp, axes):
@pmp("shp", shapes) @pmp("shp", shapes)
def test_identity(shp): def test_identity(shp):
a = np.random.rand(*shp)-0.5 + 1j*np.random.rand(*shp)-0.5j rng = np.random.default_rng(np.random.SeedSequence(42))
a = rng.random(shp)-0.5 + 1j*rng.random(shp)-0.5j
assert_(_l2error(ifftn(fftn(a), inorm=2), a) < 1.5e-15) assert_(_l2error(ifftn(fftn(a), inorm=2), a) < 1.5e-15)
assert_(_l2error(ifftn(fftn(a.real), inorm=2), a.real) < 1.5e-15) assert_(_l2error(ifftn(fftn(a.real), inorm=2), a.real) < 1.5e-15)
assert_(_l2error(fftn(ifftn(a.real), inorm=2), a.real) < 1.5e-15) assert_(_l2error(fftn(ifftn(a.real), inorm=2), a.real) < 1.5e-15)
...@@ -165,7 +171,8 @@ def test_identity(shp): ...@@ -165,7 +171,8 @@ def test_identity(shp):
@pmp("shp", shapes) @pmp("shp", shapes)
def test_identity_r(shp): def test_identity_r(shp):
a = np.random.rand(*shp)-0.5 rng = np.random.default_rng(np.random.SeedSequence(42))
a = rng.random(shp)-0.5
b = a.astype(np.float32) b = a.astype(np.float32)
for ax in range(a.ndim): for ax in range(a.ndim):
n = a.shape[ax] n = a.shape[ax]
...@@ -177,14 +184,16 @@ def test_identity_r(shp): ...@@ -177,14 +184,16 @@ def test_identity_r(shp):
@pmp("shp", shapes) @pmp("shp", shapes)
def test_identity_r2(shp): def test_identity_r2(shp):
a = np.random.rand(*shp)-0.5 + 1j*np.random.rand(*shp)-0.5j rng = np.random.default_rng(np.random.SeedSequence(42))
a = rng.random(shp)-0.5 + 1j*rng.random(shp)-0.5j
a = rfftn(irfftn(a)) a = rfftn(irfftn(a))
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_genuine_hartley(shp): def test_genuine_hartley(shp):
a = np.random.rand(*shp)-0.5 rng = np.random.default_rng(np.random.SeedSequence(42))
a = rng.random(shp)-0.5
v1 = pypocketfft.genuine_hartley(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
...@@ -193,14 +202,16 @@ def test_genuine_hartley(shp): ...@@ -193,14 +202,16 @@ def test_genuine_hartley(shp):
@pmp("shp", shapes) @pmp("shp", shapes)
def test_hartley_identity(shp): def test_hartley_identity(shp):
a = np.random.rand(*shp)-0.5 rng = np.random.default_rng(np.random.SeedSequence(42))
a = rng.random(shp)-0.5
v1 = pypocketfft.separable_hartley(pypocketfft.separable_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_genuine_hartley_identity(shp): def test_genuine_hartley_identity(shp):
a = np.random.rand(*shp)-0.5 rng = np.random.default_rng(np.random.SeedSequence(42))
a = rng.random(shp)-0.5
v1 = pypocketfft.genuine_hartley(pypocketfft.genuine_hartley(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()
...@@ -212,7 +223,8 @@ def test_genuine_hartley_identity(shp): ...@@ -212,7 +223,8 @@ def test_genuine_hartley_identity(shp):
@pmp("shp", shapes2D+shapes3D) @pmp("shp", shapes2D+shapes3D)
@pmp("axes", ((0,), (1,), (0, 1), (1, 0))) @pmp("axes", ((0,), (1,), (0, 1), (1, 0)))
def test_genuine_hartley_2D(shp, axes): def test_genuine_hartley_2D(shp, axes):
a = np.random.rand(*shp)-0.5 rng = np.random.default_rng(np.random.SeedSequence(42))
a = rng.random(shp)-0.5
assert_(_l2error(pypocketfft.genuine_hartley(pypocketfft.genuine_hartley( assert_(_l2error(pypocketfft.genuine_hartley(pypocketfft.genuine_hartley(
a, axes=axes), axes=axes, inorm=2), a) < 1e-15) a, axes=axes), axes=axes, inorm=2), a) < 1e-15)
...@@ -222,7 +234,8 @@ def test_genuine_hartley_2D(shp, axes): ...@@ -222,7 +234,8 @@ def test_genuine_hartley_2D(shp, axes):
@pmp("type", [1, 2, 3, 4]) @pmp("type", [1, 2, 3, 4])
@pmp("dtype", dtypes) @pmp("dtype", dtypes)
def testdcst1D(len, inorm, type, dtype): def testdcst1D(len, inorm, type, dtype):
a = (np.random.rand(len)-0.5).astype(dtype) rng = np.random.default_rng(np.random.SeedSequence(42))
a = (rng.random(len)-0.5).astype(dtype)
eps = tol[dtype] eps = tol[dtype]
itp = (0, 1, 3, 2, 4) itp = (0, 1, 3, 2, 4)
itype = itp[type] itype = itp[type]
......
...@@ -14,7 +14,8 @@ def test_GL(params): ...@@ -14,7 +14,8 @@ def test_GL(params):
job = pysharp.sharpjob_d() job = pysharp.sharpjob_d()
nalm = ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax) nalm = ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
nalm_r = nalm*2-lmax-1 nalm_r = nalm*2-lmax-1
alm_r = np.random.uniform(-1., 1., nalm_r) rng = np.random.default_rng(np.random.SeedSequence(42))
alm_r = rng.uniform(-1., 1., nalm_r)
alm = np.empty(nalm, dtype=np.complex128) alm = np.empty(nalm, dtype=np.complex128)
alm[0:lmax+1] = alm_r[0:lmax+1] alm[0:lmax+1] = alm_r[0:lmax+1]
alm[lmax+1:] = np.sqrt(0.5)*(alm_r[lmax+1::2] + 1j*alm_r[lmax+2::2]) alm[lmax+1:] = np.sqrt(0.5)*(alm_r[lmax+1::2] + 1j*alm_r[lmax+2::2])
...@@ -32,7 +33,8 @@ def test_fejer1(params): ...@@ -32,7 +33,8 @@ def test_fejer1(params):
job = pysharp.sharpjob_d() job = pysharp.sharpjob_d()
nalm = ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax) nalm = ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
nalm_r = nalm*2-lmax-1 nalm_r = nalm*2-lmax-1
alm_r = np.random.uniform(-1., 1., nalm_r) rng = np.random.default_rng(np.random.SeedSequence(42))
alm_r = rng.uniform(-1., 1., nalm_r)
alm = np.empty(nalm, dtype=np.complex128) alm = np.empty(nalm, dtype=np.complex128)
alm[0:lmax+1] = alm_r[0:lmax+1] alm[0:lmax+1] = alm_r[0:lmax+1]
alm[lmax+1:] = np.sqrt(0.5)*(alm_r[lmax+1::2] + 1j*alm_r[lmax+2::2]) alm[lmax+1:] = np.sqrt(0.5)*(alm_r[lmax+1::2] + 1j*alm_r[lmax+2::2])
...@@ -50,7 +52,8 @@ def test_dh(params): ...@@ -50,7 +52,8 @@ def test_dh(params):
job = pysharp.sharpjob_d() job = pysharp.sharpjob_d()
nalm = ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax) nalm = ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
nalm_r = nalm*2-lmax-1 nalm_r = nalm*2-lmax-1
alm_r = np.random.uniform(-1., 1., nalm_r) rng = np.random.default_rng(np.random.SeedSequence(42))
alm_r = rng.uniform(-1., 1., nalm_r)
alm = np.empty(nalm, dtype=np.complex128) alm = np.empty(nalm, dtype=np.complex128)
alm[0:lmax+1] = alm_r[0:lmax+1] alm[0:lmax+1] = alm_r[0:lmax+1]
alm[lmax+1:] = np.sqrt(0.5)*(alm_r[lmax+1::2] + 1j*alm_r[lmax+2::2]) alm[lmax+1:] = np.sqrt(0.5)*(alm_r[lmax+1::2] + 1j*alm_r[lmax+2::2])
......
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