Commit 6e056ee2 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

merge

parents 2a782e37 b7a416e3
Pipeline #76256 passed with stages
in 13 minutes and 10 seconds
[build-system]
requires = ["setuptools >= 40.6.0", "pybind11 >= 2.5.0", "numpy >= 1.17.0"]
build-backend = "setuptools.build_meta"
...@@ -3,7 +3,8 @@ import ducc_0_1.fft as duccfft ...@@ -3,7 +3,8 @@ import ducc_0_1.fft as duccfft
from time import time from time import time
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
np.random.seed(42)
rng = np.random.default_rng(42)
def _l2error(a, b): def _l2error(a, b):
...@@ -107,11 +108,11 @@ def bench_nd(ndim, nmax, nthr, ntry, tp, funcs, nrepeat, ttl="", filename="", ...@@ -107,11 +108,11 @@ def bench_nd(ndim, nmax, nthr, ntry, tp, funcs, nrepeat, ttl="", filename="",
print("{}D, type {}, max extent is {}:".format(ndim, tp, nmax)) print("{}D, type {}, max extent is {}:".format(ndim, tp, nmax))
results = [[] for i in range(len(funcs))] results = [[] for i in range(len(funcs))]
for n in range(ntry): for n in range(ntry):
shp = np.random.randint(nmax//3, nmax+1, ndim) shp = rng.integers(nmax//3, nmax+1, ndim)
if nice_sizes: if nice_sizes:
shp = np.array([duccfft.good_size(sz) for sz in shp]) shp = np.array([duccfft.good_size(sz) for sz in shp])
print(" {0:4d}/{1}: shape={2} ...".format(n, ntry, shp), end=" ", flush=True) print(" {0:4d}/{1}: shape={2} ...".format(n, ntry, shp), end=" ", flush=True)
a = (np.random.rand(*shp)-0.5 + 1j*(np.random.rand(*shp)-0.5)).astype(tp) a = (rng.random(shp)-0.5 + 1j*(rng.random(shp)-0.5)).astype(tp)
output=[] output=[]
for func, res in zip(funcs, results): for func, res in zip(funcs, results):
tmp = func(a, nrepeat, nthr) tmp = func(a, nrepeat, nthr)
......
...@@ -2,6 +2,9 @@ import numpy as np ...@@ -2,6 +2,9 @@ import numpy as np
import ducc_0_1.fft as fft import ducc_0_1.fft as fft
rng = np.random.default_rng(42)
def _l2error(a, b, axes): def _l2error(a, b, axes):
return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2))/np.log2(np.max([2,np.prod(np.take(a.shape,axes))])) return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2))/np.log2(np.max([2,np.prod(np.take(a.shape,axes))]))
...@@ -41,15 +44,15 @@ def update_err(err, name, value, shape): ...@@ -41,15 +44,15 @@ def update_err(err, name, value, shape):
def test(err): def test(err):
ndim = np.random.randint(1, 5) ndim = rng.integers(1, 5)
axlen = int((2**20)**(1./ndim)) axlen = int((2**20)**(1./ndim))
shape = np.random.randint(1, axlen, ndim) shape = rng.integers(1, axlen, ndim)
axes = np.arange(ndim) axes = np.arange(ndim)
np.random.shuffle(axes) rng.shuffle(axes)
nax = np.random.randint(1, ndim+1) nax = rng.integers(1, ndim+1)
axes = axes[:nax] axes = axes[:nax]
lastsize = shape[axes[-1]] lastsize = shape[axes[-1]]
a = np.random.rand(*shape)-0.5 + 1j*np.random.rand(*shape)-0.5j a = rng.random(shape)-0.5 + 1j*rng.random(shape)-0.5j
a_32 = a.astype(np.complex64) a_32 = a.astype(np.complex64)
b = ifftn(fftn(a, axes=axes, nthreads=nthreads), axes=axes, inorm=2, b = ifftn(fftn(a, axes=axes, nthreads=nthreads), axes=axes, inorm=2,
nthreads=nthreads) nthreads=nthreads)
......
...@@ -3,17 +3,19 @@ import math ...@@ -3,17 +3,19 @@ import math
import numpy as np import numpy as np
import ducc_0_1.healpix as ph import ducc_0_1.healpix as ph
rng = np.random.default_rng(42)
def report (name,vlen,ntry,nside,isnest,perf): def report (name,vlen,ntry,nside,isnest,perf):
print (name,": ",perf*1e-6,"MOps/s",sep="") print (name,": ",perf*1e-6,"MOps/s",sep="")
def random_ptg(vlen): def random_ptg(vlen):
res=np.empty((vlen,2),dtype=np.float64) res=np.empty((vlen,2),dtype=np.float64)
res[:,0]=np.arccos((np.random.random_sample(vlen)-0.5)*2) res[:,0]=np.arccos((rng.random(vlen)-0.5)*2)
res[:,1]=np.random.random_sample(vlen)*2*math.pi res[:,1]=rng.random(vlen)*2*math.pi
return res return res
def random_pix(nside,vlen): def random_pix(nside,vlen):
return np.random.randint(low=0,high=12*nside*nside-1,size=vlen,dtype=np.int64) return rng.integers(low=0,high=12*nside*nside-1,size=vlen,dtype=np.int64)
def dummy(vlen): def dummy(vlen):
inp=np.zeros(vlen,dtype=np.int64) inp=np.zeros(vlen,dtype=np.int64)
......
...@@ -2,10 +2,12 @@ import ducc_0_1.healpix as ph ...@@ -2,10 +2,12 @@ import ducc_0_1.healpix as ph
import numpy as np import numpy as np
import math import math
rng = np.random.default_rng(42)
def random_ptg(vlen): def random_ptg(vlen):
res = np.empty((vlen, 2), dtype=np.float64) res = np.empty((vlen, 2), dtype=np.float64)
res[:,0] = np.arccos((np.random.random_sample(vlen)-0.5)*2) res[:,0] = np.arccos((rng.random(vlen)-0.5)*2)
res[:,1] = np.random.random_sample(vlen)*2*math.pi res[:,1] = rng.random(vlen)*2*math.pi
return res return res
def check_pixangpix(vlen,ntry,nside,isnest): def check_pixangpix(vlen,ntry,nside,isnest):
...@@ -13,7 +15,7 @@ def check_pixangpix(vlen,ntry,nside,isnest): ...@@ -13,7 +15,7 @@ def check_pixangpix(vlen,ntry,nside,isnest):
cnt = 0 cnt = 0
while cnt < ntry: while cnt < ntry:
cnt += 1 cnt += 1
inp = np.random.randint(low=0, high=12*nside*nside-1, size=vlen) inp = rng.integers(low=0, high=12*nside*nside-1, size=vlen)
out = base.ang2pix(base.pix2ang(inp)) out = base.ang2pix(base.pix2ang(inp))
if not np.array_equal(inp, out): if not np.array_equal(inp, out):
raise ValueError("Test failed") raise ValueError("Test failed")
...@@ -33,7 +35,7 @@ def check_pixangvecpix(vlen, ntry, nside, isnest): ...@@ -33,7 +35,7 @@ def check_pixangvecpix(vlen, ntry, nside, isnest):
cnt = 0 cnt = 0
while cnt < ntry: while cnt < ntry:
cnt += 1 cnt += 1
inp = np.random.randint(low=0, high=12*nside*nside-1, size=vlen) inp = rng.integers(low=0, high=12*nside*nside-1, size=vlen)
out = base.vec2pix(ph.ang2vec(base.pix2ang(inp))) out = base.vec2pix(ph.ang2vec(base.pix2ang(inp)))
if not np.array_equal(inp,out): if not np.array_equal(inp,out):
raise ValueError("Test failed") raise ValueError("Test failed")
...@@ -43,7 +45,7 @@ def check_pixvecangpix(vlen, ntry, nside, isnest): ...@@ -43,7 +45,7 @@ def check_pixvecangpix(vlen, ntry, nside, isnest):
cnt = 0 cnt = 0
while cnt < ntry: while cnt < ntry:
cnt += 1 cnt += 1
inp = np.random.randint(low=0, high=12*nside*nside-1, size=vlen) inp = rng.integers(low=0, high=12*nside*nside-1, size=vlen)
out = base.ang2pix(ph.vec2ang(base.pix2vec(inp))) out = base.ang2pix(ph.vec2ang(base.pix2vec(inp)))
if not np.array_equal(inp,out): if not np.array_equal(inp,out):
raise ValueError("Test failed") raise ValueError("Test failed")
...@@ -53,7 +55,7 @@ def check_pixvecpix(vlen,ntry,nside,isnest): ...@@ -53,7 +55,7 @@ def check_pixvecpix(vlen,ntry,nside,isnest):
cnt=0 cnt=0
while (cnt<ntry): while (cnt<ntry):
cnt+=1 cnt+=1
inp=np.random.randint(low=0,high=12*nside*nside-1,size=vlen) inp=rng.integers(low=0,high=12*nside*nside-1,size=vlen)
out=base.vec2pix(base.pix2vec(inp)) out=base.vec2pix(base.pix2vec(inp))
if (np.array_equal(inp,out)==False): if (np.array_equal(inp,out)==False):
raise ValueError("Test failed") raise ValueError("Test failed")
...@@ -63,7 +65,7 @@ def check_ringnestring(vlen,ntry,nside): ...@@ -63,7 +65,7 @@ def check_ringnestring(vlen,ntry,nside):
cnt=0 cnt=0
while (cnt<ntry): while (cnt<ntry):
cnt+=1 cnt+=1
inp=np.random.randint(low=0,high=12*nside*nside-1,size=vlen) inp=rng.integers(low=0,high=12*nside*nside-1,size=vlen)
out=base.nest2ring(base.ring2nest(inp)) out=base.nest2ring(base.ring2nest(inp))
if (np.array_equal(inp,out)==False): if (np.array_equal(inp,out)==False):
raise ValueError("Test failed") raise ValueError("Test failed")
...@@ -73,7 +75,7 @@ def check_pixxyfpix(vlen,ntry,nside,isnest): ...@@ -73,7 +75,7 @@ def check_pixxyfpix(vlen,ntry,nside,isnest):
cnt=0 cnt=0
while (cnt<ntry): while (cnt<ntry):
cnt+=1 cnt+=1
inp=np.random.randint(low=0,high=12*nside*nside-1,size=vlen) inp=rng.integers(low=0,high=12*nside*nside-1,size=vlen)
out=base.xyf2pix(base.pix2xyf(inp)) out=base.xyf2pix(base.pix2xyf(inp))
if (np.array_equal(inp,out)==False): if (np.array_equal(inp,out)==False):
raise ValueError("Test failed") raise ValueError("Test failed")
......
...@@ -35,7 +35,8 @@ job = sht.sharpjob_d() ...@@ -35,7 +35,8 @@ job = sht.sharpjob_d()
# number of required a_lm coefficients # number of required a_lm coefficients
nalm = ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax) nalm = ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
# get random a_lm # get random a_lm
alm = np.random.uniform(-1., 1., nalm) + 1j*np.random.uniform(-1., 1., nalm) rng = np.random.default_rng(42)
alm = rng.uniform(-1., 1., nalm) + 1j*rng.uniform(-1., 1., nalm)
# make a_lm with m==0 real-valued # make a_lm with m==0 real-valued
alm[0:lmax+1].imag = 0. alm[0:lmax+1].imag = 0.
......
...@@ -14,7 +14,8 @@ print("Generating spherical harmonic coefficients up to {}".format(lmax)) ...@@ -14,7 +14,8 @@ print("Generating spherical harmonic coefficients up to {}".format(lmax))
# number of required a_lm coefficients # number of required a_lm coefficients
nalm = ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax) nalm = ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
# get random a_lm # get random a_lm
alm = np.random.uniform(-1., 1., nalm) + 1j*np.random.uniform(-1., 1., nalm) rng = np.random.default_rng(42)
alm = rng.uniform(-1., 1., nalm) + 1j*rng.uniform(-1., 1., nalm)
# make a_lm with m==0 real-valued # make a_lm with m==0 real-valued
alm[0:lmax+1].imag = 0. alm[0:lmax+1].imag = 0.
......
...@@ -5,15 +5,15 @@ import ducc_0_1.misc as misc ...@@ -5,15 +5,15 @@ import ducc_0_1.misc as misc
import time import time
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
np.random.seed(48) rng = np.random.default_rng(42)
def nalm(lmax, mmax): def nalm(lmax, mmax):
return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax) return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
def random_alm(lmax, mmax, ncomp): def random_alm(lmax, mmax, ncomp):
res = np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \ res = rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \
+ 1j*np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) + 1j*rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp))
# make a_lm with m==0 real-valued # make a_lm with m==0 real-valued
res[0:lmax+1,:].imag = 0. res[0:lmax+1,:].imag = 0.
return res return res
......
...@@ -2,15 +2,15 @@ import ducc_0_1.totalconvolve as totalconvolve ...@@ -2,15 +2,15 @@ import ducc_0_1.totalconvolve as totalconvolve
import numpy as np import numpy as np
import time import time
np.random.seed(48) rng = np.random.default_rng(48)
def nalm(lmax, mmax): def nalm(lmax, mmax):
return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax) return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
def random_alm(lmax, mmax, ncomp): def random_alm(lmax, mmax, ncomp):
res = np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \ res = rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \
+ 1j*np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) + 1j*rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp))
# make a_lm with m==0 real-valued # make a_lm with m==0 real-valued
res[0:lmax+1,:].imag = 0. res[0:lmax+1,:].imag = 0.
return res return res
...@@ -68,7 +68,7 @@ t0=time.time() ...@@ -68,7 +68,7 @@ t0=time.time()
nth = lmax+1 nth = lmax+1
nph = 2*lmax+1 nph = 2*lmax+1
ptg=np.random.uniform(0.,1.,3*nptg).reshape(nptg,3) ptg=rng.uniform(0.,1.,3*nptg).reshape(nptg,3)
ptg[:,0]*=np.pi ptg[:,0]*=np.pi
ptg[:,1]*=2*np.pi ptg[:,1]*=2*np.pi
ptg[:,2]*=2*np.pi ptg[:,2]*=2*np.pi
...@@ -78,7 +78,7 @@ bar=foo.interpol(ptg) ...@@ -78,7 +78,7 @@ bar=foo.interpol(ptg)
del foo del foo
print("Interpolating {} random angle triplets: {}s".format(nptg, time.time() -t0)) print("Interpolating {} random angle triplets: {}s".format(nptg, time.time() -t0))
t0=time.time() t0=time.time()
fake = np.random.uniform(0.,1., (ptg.shape[0],ncomp2)) fake = rng.uniform(0.,1., (ptg.shape[0],ncomp2))
foo2 = totalconvolve.PyInterpolator(lmax, kmax, ncomp2, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads) foo2 = totalconvolve.PyInterpolator(lmax, kmax, ncomp2, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
t0=time.time() t0=time.time()
foo2.deinterpol(ptg.reshape((-1,3)), fake) foo2.deinterpol(ptg.reshape((-1,3)), fake)
......
...@@ -24,6 +24,8 @@ using namespace mr; ...@@ -24,6 +24,8 @@ using namespace mr;
PYBIND11_MODULE(PKGNAME, m) PYBIND11_MODULE(PKGNAME, m)
{ {
m.attr("__version__") = PKGVERSION;
add_fft(m); add_fft(m);
add_sht(m); add_sht(m);
add_totalconvolve(m); add_totalconvolve(m);
......
...@@ -71,7 +71,7 @@ dtypes = [np.float32, np.float64, ...@@ -71,7 +71,7 @@ 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):
rng = np.random.default_rng(np.random.SeedSequence(42)) rng = np.random.default_rng(42)
a = rng.random(len)-0.5 + 1j*rng.random(len)-0.5j 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]
...@@ -92,7 +92,7 @@ def test1D(len, inorm, dtype): ...@@ -92,7 +92,7 @@ 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):
rng = np.random.default_rng(np.random.SeedSequence(42)) rng = np.random.default_rng(42)
a = rng.random(shp)-0.5 + 1j*rng.random(shp)-0.5j 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)
...@@ -105,7 +105,7 @@ def test_fftn(shp, nthreads, inorm): ...@@ -105,7 +105,7 @@ 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):
rng = np.random.default_rng(np.random.SeedSequence(42)) rng = np.random.default_rng(42)
a = rng.random(shp)-0.5 + 1j*rng.random(shp)-0.5j 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)
...@@ -116,7 +116,7 @@ def test_fftn2D(shp, axes, inorm): ...@@ -116,7 +116,7 @@ def test_fftn2D(shp, axes, inorm):
@pmp("shp", shapes) @pmp("shp", shapes)
def test_rfftn(shp): def test_rfftn(shp):
rng = np.random.default_rng(np.random.SeedSequence(42)) rng = np.random.default_rng(42)
a = rng.random(shp)-0.5 a = rng.random(shp)-0.5
tmp1 = rfftn(a) tmp1 = rfftn(a)
tmp2 = fftn(a) tmp2 = fftn(a)
...@@ -142,7 +142,7 @@ def test_rfftn(shp): ...@@ -142,7 +142,7 @@ 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):
rng = np.random.default_rng(np.random.SeedSequence(42)) rng = np.random.default_rng(42)
a = rng.random(shp)-0.5 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)
...@@ -157,7 +157,7 @@ def test_rfftn2D(shp, axes): ...@@ -157,7 +157,7 @@ def test_rfftn2D(shp, axes):
@pmp("shp", shapes) @pmp("shp", shapes)
def test_identity(shp): def test_identity(shp):
rng = np.random.default_rng(np.random.SeedSequence(42)) rng = np.random.default_rng(42)
a = rng.random(shp)-0.5 + 1j*rng.random(shp)-0.5j 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)
...@@ -171,7 +171,7 @@ def test_identity(shp): ...@@ -171,7 +171,7 @@ def test_identity(shp):
@pmp("shp", shapes) @pmp("shp", shapes)
def test_identity_r(shp): def test_identity_r(shp):
rng = np.random.default_rng(np.random.SeedSequence(42)) rng = np.random.default_rng(42)
a = rng.random(shp)-0.5 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):
...@@ -184,7 +184,7 @@ def test_identity_r(shp): ...@@ -184,7 +184,7 @@ def test_identity_r(shp):
@pmp("shp", shapes) @pmp("shp", shapes)
def test_identity_r2(shp): def test_identity_r2(shp):
rng = np.random.default_rng(np.random.SeedSequence(42)) rng = np.random.default_rng(42)
a = rng.random(shp)-0.5 + 1j*rng.random(shp)-0.5j 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)
...@@ -192,7 +192,7 @@ def test_identity_r2(shp): ...@@ -192,7 +192,7 @@ def test_identity_r2(shp):
@pmp("shp", shapes2D+shapes3D) @pmp("shp", shapes2D+shapes3D)
def test_genuine_hartley(shp): def test_genuine_hartley(shp):
rng = np.random.default_rng(np.random.SeedSequence(42)) rng = np.random.default_rng(42)
a = rng.random(shp)-0.5 a = rng.random(shp)-0.5
v1 = fft.genuine_hartley(a) v1 = fft.genuine_hartley(a)
v2 = fftn(a.astype(np.complex128)) v2 = fftn(a.astype(np.complex128))
...@@ -202,7 +202,7 @@ def test_genuine_hartley(shp): ...@@ -202,7 +202,7 @@ def test_genuine_hartley(shp):
@pmp("shp", shapes) @pmp("shp", shapes)
def test_hartley_identity(shp): def test_hartley_identity(shp):
rng = np.random.default_rng(np.random.SeedSequence(42)) rng = np.random.default_rng(42)
a = rng.random(shp)-0.5 a = rng.random(shp)-0.5
v1 = fft.separable_hartley(fft.separable_hartley(a))/a.size v1 = fft.separable_hartley(fft.separable_hartley(a))/a.size
assert_(_l2error(a, v1) < 1e-15) assert_(_l2error(a, v1) < 1e-15)
...@@ -210,7 +210,7 @@ def test_hartley_identity(shp): ...@@ -210,7 +210,7 @@ def test_hartley_identity(shp):
@pmp("shp", shapes) @pmp("shp", shapes)
def test_genuine_hartley_identity(shp): def test_genuine_hartley_identity(shp):
rng = np.random.default_rng(np.random.SeedSequence(42)) rng = np.random.default_rng(42)
a = rng.random(shp)-0.5 a = rng.random(shp)-0.5
v1 = fft.genuine_hartley(fft.genuine_hartley(a))/a.size v1 = fft.genuine_hartley(fft.genuine_hartley(a))/a.size
assert_(_l2error(a, v1) < 1e-15) assert_(_l2error(a, v1) < 1e-15)
...@@ -223,7 +223,7 @@ def test_genuine_hartley_identity(shp): ...@@ -223,7 +223,7 @@ 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):
rng = np.random.default_rng(np.random.SeedSequence(42)) rng = np.random.default_rng(42)
a = rng.random(shp)-0.5 a = rng.random(shp)-0.5
assert_(_l2error(fft.genuine_hartley(fft.genuine_hartley( assert_(_l2error(fft.genuine_hartley(fft.genuine_hartley(
a, axes=axes), axes=axes, inorm=2), a) < 1e-15) a, axes=axes), axes=axes, inorm=2), a) < 1e-15)
...@@ -234,7 +234,7 @@ def test_genuine_hartley_2D(shp, axes): ...@@ -234,7 +234,7 @@ 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):
rng = np.random.default_rng(np.random.SeedSequence(42)) rng = np.random.default_rng(42)
a = (rng.random(len)-0.5).astype(dtype) 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)
......
...@@ -23,50 +23,56 @@ nside_ring = list2fixture(pow2+nonpow2) ...@@ -23,50 +23,56 @@ nside_ring = list2fixture(pow2+nonpow2)
vlen = list2fixture([1, 10, 100, 1000, 10000]) vlen = list2fixture([1, 10, 100, 1000, 10000])
def random_ptg(vlen): def random_ptg(rng, vlen):
res = np.empty((vlen, 2), dtype=np.float64) res = np.empty((vlen, 2), dtype=np.float64)
res[:, 0] = np.arccos((np.random.random_sample(vlen)-0.5)*2) res[:, 0] = np.arccos((rng.random(vlen)-0.5)*2)
# res[:, 0] = math.pi*np.random.random_sample(vlen) # res[:, 0] = math.pi*rng.random(vlen)
res[:, 1] = np.random.random_sample(vlen)*2*math.pi res[:, 1] = rng.random(vlen)*2*math.pi
return res return res
def test_pixangpix_nest(vlen, nside_nest): def test_pixangpix_nest(vlen, nside_nest):
base = ph.Healpix_Base(nside_nest, "NEST") base = ph.Healpix_Base(nside_nest, "NEST")
inp = np.random.randint(low=0, high=12*nside_nest*nside_nest-1, size=vlen) rng = np.random.default_rng(42)
inp = rng.integers(low=0, high=12*nside_nest*nside_nest-1, size=vlen)
out = base.ang2pix(base.pix2ang(inp)) out = base.ang2pix(base.pix2ang(inp))
assert_equal(inp, out) assert_equal(inp, out)
def test_pixangpix_ring(vlen, nside_ring): def test_pixangpix_ring(vlen, nside_ring):
base = ph.Healpix_Base(nside_ring, "RING") base = ph.Healpix_Base(nside_ring, "RING")
inp = np.random.randint(low=0, high=12*nside_ring*nside_ring-1, size=vlen) rng = np.random.default_rng(42)
inp = rng.integers(low=0, high=12*nside_ring*nside_ring-1, size=vlen)
out = base