use new numpy random interface everywhere

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)