Merge branch 'ducc_0_1' into tweak_setup

python/demos/fft_bench.py

@@ -3,7 +3,8 @@ import ducc_0_1.fft as duccfft
 from time import time
 import matplotlib.pyplot as plt
 
-np.random.seed(42)
+
+rng = np.random.default_rng(42)
 
 
 def _l2error(a, b):
@@ -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))
     results = [[] for i in range(len(funcs))]
     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:
             shp = np.array([duccfft.good_size(sz) for sz in shp])
         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=[]
         for func, res in zip(funcs, results):
             tmp = func(a, nrepeat, nthr)

python/demos/fft_stress.py

@@ -2,6 +2,9 @@ import numpy as np
 import ducc_0_1.fft as fft
 
 
+rng = np.random.default_rng(42)
+
+
 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))]))
 
@@ -41,15 +44,15 @@ def update_err(err, name, value, shape):
 
 
 def test(err):
-    ndim = np.random.randint(1, 5)
+    ndim = rng.integers(1, 5)
     axlen = int((2**20)**(1./ndim))
-    shape = np.random.randint(1, axlen, ndim)
+    shape = rng.integers(1, axlen, ndim)
     axes = np.arange(ndim)
-    np.random.shuffle(axes)
-    nax = np.random.randint(1, ndim+1)
+    rng.shuffle(axes)
+    nax = rng.integers(1, ndim+1)
     axes = axes[:nax]
     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)
     b = ifftn(fftn(a, axes=axes, nthreads=nthreads), axes=axes, inorm=2,
               nthreads=nthreads)

python/demos/healpix_perftest.py

@@ -3,17 +3,19 @@ import math
 import numpy as np
 import ducc_0_1.healpix as ph
 
+rng = np.random.default_rng(42)
+
 def report (name,vlen,ntry,nside,isnest,perf):
   print (name,": ",perf*1e-6,"MOps/s",sep="")
 
 def random_ptg(vlen):
   res=np.empty((vlen,2),dtype=np.float64)
-  res[:,0]=np.arccos((np.random.random_sample(vlen)-0.5)*2)
-  res[:,1]=np.random.random_sample(vlen)*2*math.pi
+  res[:,0]=np.arccos((rng.random(vlen)-0.5)*2)
+  res[:,1]=rng.random(vlen)*2*math.pi
   return res
 
 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):
   inp=np.zeros(vlen,dtype=np.int64)

python/demos/healpix_test.py

@@ -2,10 +2,12 @@ import ducc_0_1.healpix as ph
 import numpy as np
 import math
 
+rng = np.random.default_rng(42)
+
 def random_ptg(vlen):
   res = np.empty((vlen, 2), dtype=np.float64)
-  res[:,0] = np.arccos((np.random.random_sample(vlen)-0.5)*2)
-  res[:,1] = np.random.random_sample(vlen)*2*math.pi
+  res[:,0] = np.arccos((rng.random(vlen)-0.5)*2)
+  res[:,1] = rng.random(vlen)*2*math.pi
   return res
 
 def check_pixangpix(vlen,ntry,nside,isnest):
@@ -13,7 +15,7 @@ def check_pixangpix(vlen,ntry,nside,isnest):
   cnt = 0
   while cnt < ntry:
     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))
     if not np.array_equal(inp, out):
       raise ValueError("Test failed")
@@ -33,7 +35,7 @@ def check_pixangvecpix(vlen, ntry, nside, isnest):
   cnt = 0
   while cnt < ntry:
     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)))
     if not np.array_equal(inp,out):
       raise ValueError("Test failed")
@@ -43,7 +45,7 @@ def check_pixvecangpix(vlen, ntry, nside, isnest):
   cnt = 0
   while cnt < ntry:
     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)))
     if not np.array_equal(inp,out):
       raise ValueError("Test failed")
@@ -53,7 +55,7 @@ def check_pixvecpix(vlen,ntry,nside,isnest):
   cnt=0
   while (cnt<ntry):
     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))
     if (np.array_equal(inp,out)==False):
       raise ValueError("Test failed")
@@ -63,7 +65,7 @@ def check_ringnestring(vlen,ntry,nside):
   cnt=0
   while (cnt<ntry):
     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))
     if (np.array_equal(inp,out)==False):
       raise ValueError("Test failed")
@@ -73,7 +75,7 @@ def check_pixxyfpix(vlen,ntry,nside,isnest):
   cnt=0
   while (cnt<ntry):
     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))
     if (np.array_equal(inp,out)==False):
       raise ValueError("Test failed")

python/demos/sht_demo.py

@@ -35,7 +35,8 @@ job = sht.sharpjob_d()
 # number of required a_lm coefficients
 nalm = ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
 # 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
 alm[0:lmax+1].imag = 0.
 

python/demos/sht_upsample_demo.py

@@ -14,7 +14,8 @@ print("Generating spherical harmonic coefficients up to {}".format(lmax))
 # number of required a_lm coefficients
 nalm = ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
 # 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
 alm[0:lmax+1].imag = 0.
 

python/demos/totalconvolve_accuracy.py

@@ -5,15 +5,15 @@ import ducc_0_1.misc as misc
 import time
 import matplotlib.pyplot as plt
 
-np.random.seed(48)
+rng = np.random.default_rng(42)
 
 def nalm(lmax, mmax):
     return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
 
 
 def random_alm(lmax, mmax, ncomp):
-    res = np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \
-     + 1j*np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp))
+    res = rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \
+     + 1j*rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp))
     # make a_lm with m==0 real-valued
     res[0:lmax+1,:].imag = 0.
     return res

python/demos/totalconvolve_demo.py

@@ -2,15 +2,15 @@ import ducc_0_1.totalconvolve as totalconvolve
 import numpy as np
 import time
 
-np.random.seed(48)
+rng = np.random.default_rng(48)
 
 def nalm(lmax, mmax):
     return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
 
 
 def random_alm(lmax, mmax, ncomp):
-    res = np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \
-     + 1j*np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp))
+    res = rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \
+     + 1j*rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp))
     # make a_lm with m==0 real-valued
     res[0:lmax+1,:].imag = 0.
     return res
@@ -68,7 +68,7 @@ t0=time.time()
 nth = 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[:,1]*=2*np.pi
 ptg[:,2]*=2*np.pi
@@ -78,7 +78,7 @@ bar=foo.interpol(ptg)
 del foo
 print("Interpolating {} random angle triplets: {}s".format(nptg, time.time() -t0))
 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)
 t0=time.time()
 foo2.deinterpol(ptg.reshape((-1,3)), fake)

python/test/test_fft.py

@@ -71,7 +71,7 @@ dtypes = [np.float32, np.float64,
 @pmp("inorm", [0, 1, 2])
 @pmp("dtype", dtypes)
 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 = a.astype(ctype[dtype])
     eps = tol[dtype]
@@ -92,7 +92,7 @@ def test1D(len, inorm, dtype):
 @pmp("nthreads", (0, 1, 2))
 @pmp("inorm", [0, 1, 2])
 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
     assert_(_l2error(a, ifftn(fftn(a, nthreads=nthreads, inorm=inorm),
                               nthreads=nthreads, inorm=2-inorm)) < 1e-15)
@@ -105,7 +105,7 @@ def test_fftn(shp, nthreads, inorm):
 @pmp("axes", ((0,), (1,), (0, 1), (1, 0)))
 @pmp("inorm", [0, 1, 2])
 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
     assert_(_l2error(a, ifftn(fftn(a, axes=axes, inorm=inorm),
                               axes=axes, inorm=2-inorm)) < 1e-15)
@@ -116,7 +116,7 @@ def test_fftn2D(shp, axes, inorm):
 
 @pmp("shp", shapes)
 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
     tmp1 = rfftn(a)
     tmp2 = fftn(a)
@@ -142,7 +142,7 @@ def test_rfftn(shp):
 @pmp("shp", shapes2D)
 @pmp("axes", ((0,), (1,), (0, 1), (1, 0)))
 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
     tmp1 = rfftn(a,axes=axes)
     tmp2 = fftn(a,axes=axes)
@@ -157,7 +157,7 @@ def test_rfftn2D(shp, axes):
 
 @pmp("shp", shapes)
 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
     assert_(_l2error(ifftn(fftn(a), inorm=2), a) < 1.5e-15)
     assert_(_l2error(ifftn(fftn(a.real), inorm=2), a.real) < 1.5e-15)
@@ -171,7 +171,7 @@ def test_identity(shp):
 
 @pmp("shp", shapes)
 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
     b = a.astype(np.float32)
     for ax in range(a.ndim):
@@ -184,7 +184,7 @@ def test_identity_r(shp):
 
 @pmp("shp", shapes)
 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 = rfftn(irfftn(a))
     assert_(_l2error(rfftn(irfftn(a), inorm=2), a) < 1e-15)
@@ -192,7 +192,7 @@ def test_identity_r2(shp):
 
 @pmp("shp", shapes2D+shapes3D)
 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
     v1 = fft.genuine_hartley(a)
     v2 = fftn(a.astype(np.complex128))
@@ -202,7 +202,7 @@ def test_genuine_hartley(shp):
 
 @pmp("shp", shapes)
 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
     v1 = fft.separable_hartley(fft.separable_hartley(a))/a.size
     assert_(_l2error(a, v1) < 1e-15)
@@ -210,7 +210,7 @@ def test_hartley_identity(shp):
 
 @pmp("shp", shapes)
 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
     v1 = fft.genuine_hartley(fft.genuine_hartley(a))/a.size
     assert_(_l2error(a, v1) < 1e-15)
@@ -223,7 +223,7 @@ def test_genuine_hartley_identity(shp):
 @pmp("shp", shapes2D+shapes3D)
 @pmp("axes", ((0,), (1,), (0, 1), (1, 0)))
 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
     assert_(_l2error(fft.genuine_hartley(fft.genuine_hartley(
         a, axes=axes), axes=axes, inorm=2), a) < 1e-15)
@@ -234,7 +234,7 @@ def test_genuine_hartley_2D(shp, axes):
 @pmp("type", [1, 2, 3, 4])
 @pmp("dtype", dtypes)
 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)
     eps = tol[dtype]
     itp = (0, 1, 3, 2, 4)

python/test/test_healpix.py

@@ -23,50 +23,56 @@ nside_ring = list2fixture(pow2+nonpow2)
 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[:, 0] = np.arccos((np.random.random_sample(vlen)-0.5)*2)
-#    res[:, 0] = math.pi*np.random.random_sample(vlen)
-    res[:, 1] = np.random.random_sample(vlen)*2*math.pi
+    res[:, 0] = np.arccos((rng.random(vlen)-0.5)*2)
+#    res[:, 0] = math.pi*rng.random(vlen)
+    res[:, 1] = rng.random(vlen)*2*math.pi
     return res
 
 
 def test_pixangpix_nest(vlen, nside_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))
     assert_equal(inp, out)
 
 
 def test_pixangpix_ring(vlen, nside_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))
     assert_equal(inp, out)
 
 
 def test_vecpixvec_nest(vlen, nside_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))
     assert_equal(np.all(ph.v_angle(inp, out) < base.max_pixrad()), True)
 
 
 def test_vecpixvec_ring(vlen, nside_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))
     assert_equal(np.all(ph.v_angle(inp, out) < base.max_pixrad()), True)
 
 
 def test_ringnestring(vlen, nside_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))
     assert_equal(np.all(out == inp), True)
 
 
 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))
     assert_equal(np.all(np.abs(out-inp) < 1e-14), True)

python/test/test_totalconvolve.py

@@ -23,9 +23,9 @@ def nalm(lmax, mmax):
     return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
 
 
-def random_alm(lmax, mmax, ncomp):
-    res = np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \
-     + 1j*np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp))
+def random_alm(rng, lmax, mmax, ncomp):
+    res = rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \
+     + 1j*rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp))
     # make a_lm with m==0 real-valued
     res[0:lmax+1,:].imag = 0.
     return res
@@ -57,16 +57,17 @@ def myalmdot(a1,a2,lmax,mmax,spin):
 @pmp("separate", [True, False])
 def test_against_convolution(lkmax, ncomp, separate):
     lmax, kmax = lkmax
-    slm = random_alm(lmax, lmax, ncomp)
-    blm = random_alm(lmax, kmax, ncomp)
+    rng = np.random.default_rng(42)
+    slm = random_alm(rng, lmax, lmax, ncomp)
+    blm = random_alm(rng, lmax, kmax, ncomp)
 
     inter = totalconvolve.PyInterpolator(slm, blm, separate, lmax, kmax,
                                          epsilon=1e-8, nthreads=2)
     nptg = 50
     ptg = np.zeros((nptg,3))
-    ptg[:,0] = np.random.uniform(0, np.pi, nptg)
-    ptg[:,1] = np.random.uniform(0, 2*np.pi, nptg)
-    ptg[:,2] = np.random.uniform(-np.pi, np.pi, nptg)
+    ptg[:,0] = rng.uniform(0, np.pi, nptg)
+    ptg[:,1] = rng.uniform(0, 2*np.pi, nptg)
+    ptg[:,2] = rng.uniform(-np.pi, np.pi, nptg)
 
     res1 = inter.interpol(ptg)
 
@@ -87,17 +88,18 @@ def test_against_convolution(lkmax, ncomp, separate):
 @pmp("separate", [True, False])
 def test_adjointness(lkmax, ncomp, separate):
     lmax, kmax = lkmax
-    slm = random_alm(lmax, lmax, ncomp)
-    blm = random_alm(lmax, kmax, ncomp)
+    rng = np.random.default_rng(42)
+    slm = random_alm(rng, lmax, lmax, ncomp)
+    blm = random_alm(rng, lmax, kmax, ncomp)
     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[:,1]*=2*np.pi
     ptg[:,2]*=2*np.pi
     foo = totalconvolve.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=1e-6, nthreads=2)
     inter1=foo.interpol(ptg)
     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.deinterpol(ptg.reshape((-1,3)), fake)
     bla=foo2.getSlm(blm)

python/test/test_wgridder.py

@@ -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):
     if singleprec and epsilon < 5e-5:
         return
-    np.random.seed(42)
+    rng = np.random.default_rng(42)
     pixsizex = np.pi/180/60/nxdirty*0.2398
     pixsizey = np.pi/180/60/nxdirty
     speedoflight, f0 = 299792458., 1e9
     freq = f0 + np.arange(nchan)*(f0/nchan)
-    uvw = (np.random.rand(nrow, 3)-0.5)/(pixsizey*f0/speedoflight)
-    ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5)
-    wgt = np.random.rand(nrow, nchan) if use_wgt else None
-    dirty = np.random.rand(nxdirty, nydirty)-0.5
+    uvw = (rng.random((nrow, 3))-0.5)/(pixsizey*f0/speedoflight)
+    ms = rng.random((nrow, nchan))-0.5 + 1j*(rng.random((nrow, nchan))-0.5)
+    wgt = rng.random((nrow, nchan)) if use_wgt else None
+    dirty = rng.random((nxdirty, nydirty))-0.5
     nu, nv = int(nxdirty*ofactor)+1, int(nydirty*ofactor)+1
     if nu&1:
         nu+=1
@@ -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):
     if singleprec and epsilon < 5e-5:
         return
-    np.random.seed(40)
+    rng = np.random.default_rng(42)
     pixsizex = fov*np.pi/180/nxdirty
     pixsizey = fov*np.pi/180/nydirty*1.1
     speedoflight, f0 = 299792458., 1e9
     freq = f0 + np.arange(nchan)*(f0/nchan)
-    uvw = (np.random.rand(nrow, 3)-0.5)/(pixsizex*f0/speedoflight)
-    ms = np.random.rand(nrow, nchan)-0.5 + 1j*(np.random.rand(nrow, nchan)-0.5)
-    wgt = np.random.rand(nrow, nchan) if use_wgt else None
+    uvw = (rng.random((nrow, 3))-0.5)/(pixsizex*f0/speedoflight)
+    ms = rng.random((nrow, nchan))-0.5 + 1j*(rng.random((nrow, nchan))-0.5)
+    wgt = rng.random((nrow, nchan)) if use_wgt else None
     nu, nv = int(nxdirty*ofactor)+1, int(nydirty*ofactor)+1
     if nu&1:
         nu+=1