Commit bf32424e authored by Martin Reinecke's avatar Martin Reinecke

adjust tests to new numpy RNG interface

parent bd30feab
...@@ -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.ang2pix(base.pix2ang(inp)) out = base.ang2pix(base.pix2ang(inp))
assert_equal(inp, out) assert_equal(inp, out)
def test_vecpixvec_nest(vlen, nside_nest): def test_vecpixvec_nest(vlen, nside_nest):
base = ph.Healpix_Base(nside_nest, "NEST") base = ph.Healpix_Base(nside_nest, "NEST")
inp = ph.ang2vec(random_ptg(vlen)) rng = np.random.default_rng(42)
inp = ph.ang2vec(random_ptg(rng, vlen))
out = base.pix2vec(base.vec2pix(inp)) out = base.pix2vec(base.vec2pix(inp))
assert_equal(np.all(ph.v_angle(inp, out) < base.max_pixrad()), True) assert_equal(np.all(ph.v_angle(inp, out) < base.max_pixrad()), True)
def test_vecpixvec_ring(vlen, nside_ring): def test_vecpixvec_ring(vlen, nside_ring):
base = ph.Healpix_Base(nside_ring, "RING") base = ph.Healpix_Base(nside_ring, "RING")
inp = ph.ang2vec(random_ptg(vlen)) rng = np.random.default_rng(42)
inp = ph.ang2vec(random_ptg(rng, vlen))
out = base.pix2vec(base.vec2pix(inp)) out = base.pix2vec(base.vec2pix(inp))
assert_equal(np.all(ph.v_angle(inp, out) < base.max_pixrad()), True) assert_equal(np.all(ph.v_angle(inp, out) < base.max_pixrad()), True)
def test_ringnestring(vlen, nside_nest): def test_ringnestring(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.ring2nest(base.nest2ring(inp)) out = base.ring2nest(base.nest2ring(inp))
assert_equal(np.all(out == inp), True) assert_equal(np.all(out == inp), True)
def test_vecangvec(vlen): def test_vecangvec(vlen):
inp = random_ptg(vlen) rng = np.random.default_rng(42)
inp = random_ptg(rng, vlen)
out = ph.vec2ang(ph.ang2vec(inp)) out = ph.vec2ang(ph.ang2vec(inp))
assert_equal(np.all(np.abs(out-inp) < 1e-14), True) assert_equal(np.all(np.abs(out-inp) < 1e-14), True)
...@@ -23,9 +23,9 @@ def nalm(lmax, mmax): ...@@ -23,9 +23,9 @@ 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(rng, 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
...@@ -57,16 +57,17 @@ def myalmdot(a1,a2,lmax,mmax,spin): ...@@ -57,16 +57,17 @@ def myalmdot(a1,a2,lmax,mmax,spin):
@pmp("separate", [True, False]) @pmp("separate", [True, False])
def test_against_convolution(lkmax, ncomp, separate): def test_against_convolution(lkmax, ncomp, separate):
lmax, kmax = lkmax lmax, kmax = lkmax
slm = random_alm(lmax, lmax, ncomp) rng = np.random.default_rng(42)
blm = random_alm(lmax, kmax, ncomp) slm = random_alm(rng, lmax, lmax, ncomp)
blm = random_alm(rng, lmax, kmax, ncomp)
inter = totalconvolve.PyInterpolator(slm, blm, separate, lmax, kmax, inter = totalconvolve.PyInterpolator(slm, blm, separate, lmax, kmax,
epsilon=1e-8, nthreads=2) epsilon=1e-8, nthreads=2)
nptg = 50 nptg = 50
ptg = np.zeros((nptg,3)) ptg = np.zeros((nptg,3))
ptg[:,0] = np.random.uniform(0, np.pi, nptg) ptg[:,0] = rng.uniform(0, np.pi, nptg)
ptg[:,1] = np.random.uniform(0, 2*np.pi, nptg) ptg[:,1] = rng.uniform(0, 2*np.pi, nptg)
ptg[:,2] = np.random.uniform(-np.pi, np.pi, nptg) ptg[:,2] = rng.uniform(-np.pi, np.pi, nptg)
res1 = inter.interpol(ptg) res1 = inter.interpol(ptg)
...@@ -87,17 +88,18 @@ def test_against_convolution(lkmax, ncomp, separate): ...@@ -87,17 +88,18 @@ def test_against_convolution(lkmax, ncomp, separate):
@pmp("separate", [True, False]) @pmp("separate", [True, False])
def test_adjointness(lkmax, ncomp, separate): def test_adjointness(lkmax, ncomp, separate):
lmax, kmax = lkmax lmax, kmax = lkmax
slm = random_alm(lmax, lmax, ncomp) rng = np.random.default_rng(42)
blm = random_alm(lmax, kmax, ncomp) slm = random_alm(rng, lmax, lmax, ncomp)
blm = random_alm(rng, lmax, kmax, ncomp)
nptg=100000 nptg=100000
ptg=np.random.uniform(0.,1.,nptg*3).reshape(nptg,3) ptg=rng.uniform(0.,1.,nptg*3).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
foo = totalconvolve.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=1e-6, nthreads=2) foo = totalconvolve.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=1e-6, nthreads=2)
inter1=foo.interpol(ptg) inter1=foo.interpol(ptg)
ncomp2 = inter1.shape[1] ncomp2 = inter1.shape[1]
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=1e-6, nthreads=2) foo2 = totalconvolve.PyInterpolator(lmax, kmax, ncomp2, epsilon=1e-6, nthreads=2)
foo2.deinterpol(ptg.reshape((-1,3)), fake) foo2.deinterpol(ptg.reshape((-1,3)), fake)
bla=foo2.getSlm(blm) bla=foo2.getSlm(blm)
......
...@@ -50,15 +50,15 @@ def explicit_gridder(uvw, freq, ms, wgt, nxdirty, nydirty, xpixsize, ypixsize, ...@@ -50,15 +50,15 @@ def explicit_gridder(uvw, freq, ms, wgt, nxdirty, nydirty, xpixsize, ypixsize,
def test_adjointness_ms2dirty(nxdirty, nydirty, ofactor, nrow, nchan, epsilon, singleprec, wstacking, use_wgt, nthreads): def test_adjointness_ms2dirty(nxdirty, nydirty, ofactor, nrow, nchan, epsilon, singleprec, wstacking, use_wgt, nthreads):
if singleprec and epsilon < 5e-5: if singleprec and epsilon < 5e-5:
return return
np.random.seed(42) rng = np.random.default_rng(42)
pixsizex = np.pi/180/60/nxdirty*0.2398 pixsizex = np.pi/180/60/nxdirty*0.2398
pixsizey = np.pi/180/60/nxdirty pixsizey = np.pi/180/60/nxdirty
speedoflight, f0 = 299792458., 1e9 speedoflight, f0 = 299792458., 1e9
freq = f0 + np.arange(nchan)*(f0/nchan) freq = f0 + np.arange(nchan)*(f0/nchan)
uvw = (np.random.rand(nrow, 3)-0.5)/(pixsizey*f0/speedoflight) uvw = (rng.random((nrow, 3))-0.5)/(pixsizey*f0/speedoflight)
ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5) ms = rng.random((nrow, nchan))-0.5 + 1j*(rng.random((nrow, nchan))-0.5)
wgt = np.random.rand(nrow, nchan) if use_wgt else None wgt = rng.random((nrow, nchan)) if use_wgt else None
dirty = np.random.rand(nxdirty, nydirty)-0.5 dirty = rng.random((nxdirty, nydirty))-0.5
nu, nv = int(nxdirty*ofactor)+1, int(nydirty*ofactor)+1 nu, nv = int(nxdirty*ofactor)+1, int(nydirty*ofactor)+1
if nu&1: if nu&1:
nu+=1 nu+=1
...@@ -91,14 +91,14 @@ def test_adjointness_ms2dirty(nxdirty, nydirty, ofactor, nrow, nchan, epsilon, s ...@@ -91,14 +91,14 @@ def test_adjointness_ms2dirty(nxdirty, nydirty, ofactor, nrow, nchan, epsilon, s
def test_ms2dirty_against_wdft2(nxdirty, nydirty, ofactor, nrow, nchan, epsilon, singleprec, wstacking, use_wgt, fov, nthreads): def test_ms2dirty_against_wdft2(nxdirty, nydirty, ofactor, nrow, nchan, epsilon, singleprec, wstacking, use_wgt, fov, nthreads):
if singleprec and epsilon < 5e-5: if singleprec and epsilon < 5e-5:
return return
np.random.seed(40) rng = np.random.default_rng(42)
pixsizex = fov*np.pi/180/nxdirty pixsizex = fov*np.pi/180/nxdirty
pixsizey = fov*np.pi/180/nydirty*1.1 pixsizey = fov*np.pi/180/nydirty*1.1
speedoflight, f0 = 299792458., 1e9 speedoflight, f0 = 299792458., 1e9
freq = f0 + np.arange(nchan)*(f0/nchan) freq = f0 + np.arange(nchan)*(f0/nchan)
uvw = (np.random.rand(nrow, 3)-0.5)/(pixsizex*f0/speedoflight) uvw = (rng.random((nrow, 3))-0.5)/(pixsizex*f0/speedoflight)
ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5) ms = rng.random((nrow, nchan))-0.5 + 1j*(rng.random((nrow, nchan))-0.5)
wgt = np.random.rand(nrow, nchan) if use_wgt else None wgt = rng.random((nrow, nchan)) if use_wgt else None
nu, nv = int(nxdirty*ofactor)+1, int(nydirty*ofactor)+1 nu, nv = int(nxdirty*ofactor)+1, int(nydirty*ofactor)+1
if nu&1: if nu&1:
nu+=1 nu+=1
......
Markdown is supported
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